Activated and renewed HBC clustering and comparison with scATAC-seq
Koen Van den Berge
7/22/2020
here::set_here()## File .here already exists in /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/scriptssuppressPackageStartupMessages({
library(slingshot)
library(SingleCellExperiment)
library(cowplot)
library(rgl)
library(clusterExperiment)
library(RColorBrewer)
library(ggplot2)
library(pheatmap)
library(wesanderson)
library(gridExtra)
library(irlba)
library(tidyverse)
})Import data
sds <- readRDS("../data/finalTrajectory/sling.rds")
counts <- readRDS("../data/finalTrajectory/counts_noResp_noMV.rds")
# counts <- round(counts)
sce <- readRDS("../data/finalTrajectory/sce_tradeSeq20200904.rds")
load("../data/ALL_TF.Rda")
clDatta <- readRDS("../data/finalTrajectory/dattaCl_noResp_noMV.rds")
timePoint <- readRDS("../data/finalTrajectory/timePoint_noResp_noMV.rds")Select activated and renewed HBC clusters
I will select
- activated HBCs as cells sampled in the 24h timepoint that are in the starting cluster
- renewed (‘resting’) HBCs as cells sampled in the 14D timepoint that are in the HBC cluster
cols <- c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00", "#FFFF33", "#A65628", "#F781BF", "#999999")
view3d( theta = 10, phi = 10)
# Datta labels
plot3d(reducedDim(sds), aspect = 'iso', col=cols[clDatta], alpha=.6, xlab="UMAP1", ylab="UMAP2", zlab="UMAP3", box=FALSE, top=TRUE)
plot3d(sds, add=TRUE, lwd=3, col="black")
# timepoints
plot3d(reducedDim(sds), aspect = 'iso', col=cols[timePoint], alpha=.6, xlab="UMAP1", ylab="UMAP2", zlab="UMAP3", box=FALSE, top=TRUE)
plot3d(sds, add=TRUE, lwd=3, col="black")
# cluster labels used for TI
clSds <- apply(slingClusterLabels(sds), 1, which.max)
plot3d(reducedDim(sds), aspect = 'iso', col=cols[clSds], alpha=.6, xlab="UMAP1", ylab="UMAP2", zlab="UMAP3", box=FALSE, top=TRUE)
plot3d(sds, add=TRUE, lwd=3, col="black")
actHBCId <- which(clSds == 8 & timePoint == "24HPI")
restHBCId <- which(clDatta == "HBC" & timePoint == "14DPI")
grpHBC <- rep(NA, length(c(actHBCId, restHBCId)))
grpHBC[1:length(actHBCId)] <- "activated"
grpHBC[is.na(grpHBC)] <- "resting"
grpHBC <- factor(grpHBC)
countHBC <- counts[, c(actHBCId, restHBCId)]Dimensionality reduction
library(scran)
library(uwot)
library(scater)##
## Attaching package: 'scater'## The following object is masked from 'package:clusterExperiment':
##
## plotHeatmaplogFracs <- log(sweep(countHBC+1e-5, 2, colSums(countHBC), "/"))
gv <- modelGeneVar(countHBC)
hvg <- getTopHVGs(gv, n=1000)
pc10 <- irlba::irlba(logFracs[hvg,], nv=10)
um10 <- uwot::umap(pc10$v %*% diag(pc10$d))
plot(um10, col=alpha(c("orange", "darkseagreen3")[as.numeric(grpHBC)], .75), pch=16, cex=1/2, bty='l')
legend("topright", c("activated", "resting"), pch=16, col=c("orange", "darkseagreen3"), bty="n")
Batch correction with Harmony
library(harmony)## Loading required package: Rcppharmony_embeddings <- HarmonyMatrix(logFracs, meta_data=grpHBC, do_pca=TRUE)## Harmony 1/10## Harmony 2/10## Harmony converged after 2 iterationsumHarmony <- uwot::umap(harmony_embeddings)
plot(umHarmony, col=alpha(c("orange", "darkseagreen3")[as.numeric(grpHBC)], .75), pch=16, cex=1/2, bty='l')
legend("topright", c("activated", "resting"), pch=16, col=c("orange", "darkseagreen3"), bty="n")
Integration: Processing and dimensionality reduction of each dataset separately
library(Signac)
library(Seurat)##
## Attaching package: 'Seurat'## The following object is masked from 'package:SummarizedExperiment':
##
## Assayslibrary(patchwork)##
## Attaching package: 'patchwork'## The following object is masked from 'package:cowplot':
##
## align_plotslibrary(ArchR)## Loading required package: data.table##
## Attaching package: 'data.table'## The following objects are masked from 'package:dplyr':
##
## between, first, last## The following object is masked from 'package:purrr':
##
## transpose## The following object is masked from 'package:SummarizedExperiment':
##
## shift## The following object is masked from 'package:GenomicRanges':
##
## shift## The following object is masked from 'package:IRanges':
##
## shift## The following objects are masked from 'package:S4Vectors':
##
## first, second## Loading required package: rhdf5## Loading required package: magrittr##
## Attaching package: 'magrittr'## The following object is masked from 'package:purrr':
##
## set_names## The following object is masked from 'package:tidyr':
##
## extractset.seed(1234)
addArchRGenome("mm10")## Setting default genome to Mm10.# scATAC-seq importing and preprocessing
# from 20200610_scATAC_injury_archr_samples_pairwiseHBC_v3.Rmd
archProj <- loadArchRProject("../data/scATAC/archR_samples_pairwise_v2")## Successfully loaded ArchRProject!##
## / |
## / \
## . / |.
## \\\ / |.
## \\\ / `|.
## \\\ / |.
## \ / |\
## \\#####\ / ||
## ==###########> / ||
## \\##==......\ / ||
## ______ = =|__ /__ || \\\
## ,--' ,----`-,__ ___/' --,-`-===================##========>
## \ ' ##_______ _____ ,--,__,=##,__ ///
## , __== ___,-,__,--'#' ===' `-' | ##,-/
## -,____,---' \\####\\________________,--\\_##,/
## ___ .______ ______ __ __ .______
## / \ | _ \ / || | | | | _ \
## / ^ \ | |_) | | ,----'| |__| | | |_) |
## / /_\ \ | / | | | __ | | /
## / _____ \ | |\ \\___ | `----.| | | | | |\ \\___.
## /__/ \__\ | _| `._____| \______||__| |__| | _| `._____|
## clHBCMerged <- archProj$clHBCMerged
plotEmbedding(ArchRProj = archProj, colorBy = "cellColData",
name = "clHBCMerged", embedding = "UMAP_cor")## ArchR logging to : ArchRLogs/ArchR-plotEmbedding-13a0937bd7f74-Date-2021-04-02_Time-15-09-58.log
## If there is an issue, please report to github with logFile!## Getting UMAP Embedding## ColorBy = cellColData## Plotting Embedding## 1
## ArchR logging successful to : ArchRLogs/ArchR-plotEmbedding-13a0937bd7f74-Date-2021-04-02_Time-15-09-58.log
archrUmap <- ArchR::getEmbedding(archProj, embedding = "UMAP_cor")
colnames(archrUmap) <- c("UMAP1", "UMAP2")
plotEmbedding(ArchRProj = archProj, colorBy = "cellColData",
name = "treat", embedding = "UMAP_cor")## ArchR logging to : ArchRLogs/ArchR-plotEmbedding-13a0979fab542-Date-2021-04-02_Time-15-10-00.log
## If there is an issue, please report to github with logFile!
## Getting UMAP Embedding
## ColorBy = cellColData
## Plotting Embedding
## 1
## ArchR logging successful to : ArchRLogs/ArchR-plotEmbedding-13a0979fab542-Date-2021-04-02_Time-15-10-00.log
clHBCMerged2 <- as.character(getCellColData(archProj, "clHBCMerged")$clHBCMerged)
clHBCMerged2[clHBCMerged2 == "C1_Inj"] <- "Activated"
clHBCMerged2[clHBCMerged2 == "C2_Hybrid"] <- "Hybrid"
clHBCMerged2[clHBCMerged2 == "C3_UI"] <- "Resting"
archrUmap$clHBCMerged2 <- factor(clHBCMerged2)
archrUmap$treatment <- getCellColData(archProj, "treat")$treat
peaks <- archProj@peakSet
names(peaks) <- paste0("peak",1:length(peaks))
atacPeakMat <- readRDS("../data/scATAC/peakMatArchrHbc_v2.rds")
rownames(atacPeakMat) <- names(peaks)
geneActivity <- readRDS("../data/scATAC/geneMatArchrHbc.rds")
atac <- Seurat::CreateSeuratObject(counts = atacPeakMat,
assay = "ATAC",
project = "10x_ATAC")
atac[["ACTIVITY"]] <- CreateAssayObject(counts = geneActivity)
atac$tech <- "atac"
atac$treatment <- archProj@cellColData$treat
atac$clHBCMerged <- clHBCMerged
# process gene activity
DefaultAssay(atac) <- "ACTIVITY"
atac <- FindVariableFeatures(atac)
atac <- NormalizeData(atac)
atac <- ScaleData(atac)## Centering and scaling data matrix# process peak matrix
DefaultAssay(atac) <- "ATAC"
# VariableFeatures(atac) <- names(which(Matrix::rowSums(atac) > 50))
# atac <- RunLSI(atac, n = 30, scale.max = NULL)
# atac <- RunUMAP(atac, reduction = "lsi", dims = 1:30)
atac <- FindTopFeatures(atac, min.cutoff = 'q0')
# Seurat 3
atac <- RunLSI(atac, n = 30) ## Warning: RunLSI is being moved to Signac. Equivalent functionality can be
## achieved via the Signac::RunTFIDF and Signac::RunSVD functions; for more
## information on Signac, please see https://github.com/timoast/Signac
## Warning: RunLSI is being moved to Signac. Equivalent functionality can be
## achieved via the Signac::RunTFIDF and Signac::RunSVD functions; for more
## information on Signac, please see https://github.com/timoast/Signac
## Warning: RunLSI is being moved to Signac. Equivalent functionality can be
## achieved via the Signac::RunTFIDF and Signac::RunSVD functions; for more
## information on Signac, please see https://github.com/timoast/Signac## Performing TF-IDF normalization
## Running SVD on TF-IDF matrix
## Scaling cell embeddings# Seurat 4
#atac <- RunTFIDF(atac)
#atac <- RunSVD(atac, n = 30, reduction.key = "lsi_", reduction.name = "lsi")
# drop first dim as often correlated with seq depth
atac <- RunUMAP(object = atac, reduction = 'lsi', dims = 2:30)## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session## 15:11:23 UMAP embedding parameters a = 0.9922 b = 1.112
## 15:11:23 Read 4076 rows and found 29 numeric columns
## 15:11:23 Using Annoy for neighbor search, n_neighbors = 30
## 15:11:23 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 15:11:23 Writing NN index file to temp file /var/folders/15/gs2v829x4gxcdv0tlgk_5l4m0000gn/T//RtmpPH7YMA/file13a093a923e9e
## 15:11:23 Searching Annoy index using 1 thread, search_k = 3000
## 15:11:24 Annoy recall = 100%
## 15:11:26 Commencing smooth kNN distance calibration using 1 thread
## 15:11:27 Initializing from normalized Laplacian + noise
## 15:11:28 Commencing optimization for 500 epochs, with 152690 positive edges
## 15:11:34 Optimization finished# scRNA-seq import
rna <- CreateSeuratObject(counts=countHBC)
rna$tech <- '10x'
rna$celltype <- grpHBC
rna <- FindVariableFeatures(rna, selection.method = "vst", nfeatures = 2000)
rna <- ScaleData(rna)## Centering and scaling data matrixrna <- RunPCA(rna, npcs=20)## PC_ 1
## Positive: Slc6a6, Rps27, Junb, Aqp3, Btg2, Fos, Fxyd3, Apoe, Krt5, Jun
## Gclm, Socs3, Ier2, Gpm6a, Ubc, Ctsl, Cyr61, Ifitm3, Laptm4a, Atf3
## Gstm2, Csrp1, Ptn, Pltp, Fosb, Dusp1, Igfbp2, Klf4, Krt14, Gas6
## Negative: Tmsb10, Cystm1, Atf5, Plaur, Adam8, Krt6a, Map1b, Tuba1a, Tubb2b, Cst6
## Ptma, Ecm1, Stmn1, Emp3, Uchl1, Upk3bl, Ceacam1, Clcf1, Tm4sf1, Gprc5a
## Cldn7, Tubb6, Lypd3, Cldn3, Lamc2, Phlda2, Ncam1, Cldn4, Mal, Sult2b1
## PC_ 2
## Positive: Igfbp2, Cald1, Tmsb10, Ogn, Tmsb4x, Igfbp4, Adam8, Krt5, Apoe, Plaur
## Fam101b, Tspan6, Abca4, Sparc, Tm4sf1, Col23a1, Lamc2, Krt6a, Clcf1, Tmem108
## Slc29a1, Krt14, Tubb6, Ngf, Cst6, Oaf, Tubb2b, Ecm1, Hmga2, Car12
## Negative: Prdx6, Adh7, Gsta3, Sftpd, Cyp2f2, Ly6a, Tst, Slc16a11, Ifitm1, Reg3g
## Muc1, Por, Cbr2, Adh1, Lrrc26, Selenbp1, Mgst1, Gpx2, Sult1d1, Gstp1
## Cbr3, Cd14, Cyp4b1, Cxcl17, Gsto1, Serpinb11, Elf3, Ly6c1, Serpina3n, Epas1
## PC_ 3
## Positive: Adh7, Aldh3a1, Selenbp1, Adh1, Ly6a, Oat, Defb1, Fmo2, Cbr3, Gstp1
## Tmem176a, Atp1b1, Tmem176b, Efemp1, Igfbp4, Upk1b, Sparc, Ly6c1, Tgm2, Igfbp2
## Mgp, Kcnj16, Sult1d1, Gmpr, Cdo1, Krt15, Nrtn, Cyp4b1, Muc1, Bsg
## Negative: Cldn4, X1810011O10Rik, Klf4, Sfn, Zfp36, S100a10, Mast4, Atf3, Tagln2, Fosb
## Avpi1, Arc, Nr4a1, Cldn7, Nfat5, Egr1, Klf6, Cdkn1a, Phlda1, Lmna
## Crym, Neat1, Anxa1, Krt18, S100a6, Acsm4, Lypd3, Papss2, Ier2, Krt23
## PC_ 4
## Positive: mt.Co1, mt.Co3, mt.Atp6, Hpgd, Pla2r1, X1500009L16Rik, Acsm4, Rps27rt, Car12, Dapl1
## Myo6, Fth1, Papss2, Lypd2, Insl6, Cyb5r3, Scd2, Gas6, Trp63, Nmb
## Eya1, Piezo2, Limd2, Crym, Sema5a, Serpinb1b, Ftl1, X4631405K08Rik, Fndc1, Scgb1c1
## Negative: Actg1, Anxa3, Ier3, Nfkbia, F3, Tmem176b, Krt7, Sfn, Cyr61, Cldn4
## Krt14, Tacstd2, Tuba1c, Tmem176a, Fos, Bhlhe40, Junb, Tnfrsf12a, Gmpr, Cish
## Fst, Cxcl16, Plaur, Krt17, Ubc, Trib1, Oat, Atf3, Adam8, Dusp5
## PC_ 5
## Positive: Fosb, Egr1, Neat1, Irs2, Pdcd4, Nr4a1, Lmna, Klf9, Cbx3, Dst
## Ccdc3, Mast4, Sdc4, Timp3, Nfia, Kdm6b, Birc5, Maff, Spry2, Il6st
## Stmn1, Hist1h2ap, Nfat5, Sik1, Arid5a, Gm26532, Pbk, Gclc, Tpm1, Fam101b
## Negative: Aldh2, Gstm1, Prdx1, Aqp3, Ces1d, Tst, Gsn, Clu, Ephx1, Cpe
## Ifitm3, Ubc, Krt18, Prdx6, Hspb1, Galm, Fmo6, Hspa5, Ppa1, Krt5
## Krt8, Pltp, Mgst1, Ctsl, Cyp2a5, Alas1, Cryl1, Calml4, Gstm2, Actg1rna <- RunUMAP(rna, dims = 1:20, min_dist=.1)## Warning: The following arguments are not used: min_dist## 15:11:38 UMAP embedding parameters a = 0.9922 b = 1.112
## 15:11:38 Read 3064 rows and found 20 numeric columns
## 15:11:38 Using Annoy for neighbor search, n_neighbors = 30
## 15:11:38 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 15:11:38 Writing NN index file to temp file /var/folders/15/gs2v829x4gxcdv0tlgk_5l4m0000gn/T//RtmpPH7YMA/file13a09574d151c
## 15:11:38 Searching Annoy index using 1 thread, search_k = 3000
## 15:11:39 Annoy recall = 100%
## 15:11:40 Commencing smooth kNN distance calibration using 1 thread
## 15:11:42 Initializing from normalized Laplacian + noise
## 15:11:42 Commencing optimization for 500 epochs, with 122146 positive edges
## 15:11:47 Optimization finishedp1 <- DimPlot(atac, reduction = "umap", group.by="treatment", label=TRUE) +
NoLegend() +
ggtitle("scATAC-seq")## Warning: Using `as.character()` on a quosure is deprecated as of rlang 0.3.0.
## Please use `as_label()` or `as_name()` instead.
## This warning is displayed once per session.p2 <- DimPlot(rna, group.by = "celltype", label = TRUE, repel = TRUE) + NoLegend() + ggtitle("scRNA-seq")
p1 + p2
# Markers for activated/resting HBC in scRNA-seq
FeaturePlot(rna, features = "Krtdap")
FeaturePlot(rna, features = "Krt6a")
FeaturePlot(rna, features = "Krt5")
FeaturePlot(rna, features = "Krt14")
FeaturePlot(rna, features = "Trp63")
# Markers for activated/resting HBC in scATAC-seq
FeaturePlot(atac, features = "Krtdap") + xlim(c(-5,7))## Warning: Could not find Krtdap in the default search locations, found in
## ACTIVITY assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 14 rows containing missing values (geom_point).
FeaturePlot(atac, features = "Krt6a") + xlim(c(-5,7))## Warning: Could not find Krt6a in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 14 rows containing missing values (geom_point).
FeaturePlot(atac, features = "Krt5") + xlim(c(-6,6))## Warning: Could not find Krt5 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 55 rows containing missing values (geom_point).
FeaturePlot(atac, features = "Krt14") + xlim(c(-6,6))## Warning: Could not find Krt14 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 55 rows containing missing values (geom_point).
FeaturePlot(atac, features = "Trp63") + xlim(c(-6,6))## Warning: Could not find Trp63 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 55 rows containing missing values (geom_point).
Seurat reduced dimension plots
# RNA
pRNA1 <- DimPlot(rna, group.by = "celltype", label = FALSE, repel = TRUE) +
ggtitle("scRNA-seq: 24h HBC* and rHBC cells") +
scale_color_manual(labels = c("Activated", "Regenerated\n(resting)"),
values = c("steelblue", "salmon")) +
labs(color = "Cell state")
# ggsave("../plots/scRNA24hHBC.png", width=9)
# ATAC based on Seurat
pATAC1 <- DimPlot(atac, reduction = "umap", group.by="clHBCMerged", label=FALSE) +
ggtitle("scATAC-seq HBC cell states") +
xlim(c(-6, 7))
pATAC2 <- DimPlot(atac, reduction = "umap", group.by="treatment", label=FALSE, pt.size = .2) +
scale_color_manual(labels = c("Injured", "Uninjured"),
values = c("orange", "darkseagreen3")) +
ggtitle("scATAC-seq HBC treatment") +
xlim(c(-6, 7))
cowplot::plot_grid(pATAC1, pATAC2)## Warning: Removed 14 rows containing missing values (geom_point).
## Warning: Removed 14 rows containing missing values (geom_point).
# ggsave("../plots/ATACClusters_seurat.png", width=10)Data integration
Cell type transfer
# find anchor features
transfer.anchors <- FindTransferAnchors(reference = rna, query = atac,
features = VariableFeatures(object = rna),
reference.assay = "RNA",
query.assay = "ACTIVITY",
reduction = "cca")## Warning in RunCCA.Seurat(object1 = reference, object2 = query, features =
## features, : Running CCA on different assays## Running CCA## Merging objects## Finding neighborhoods## Finding anchors## Found 10751 anchors## Filtering anchors## Retained 1185 anchors## Extracting within-dataset neighbors# predict ATAC cell type based on RNA labels
celltype.predictions <- TransferData(anchorset = transfer.anchors,
refdata = rna$celltype,
weight.reduction = atac[["lsi"]])## Finding integration vectors## Finding integration vector weights## Predicting cell labels# add predicted labels to ATAC
atac <- AddMetaData(atac, metadata = celltype.predictions)
hist(atac$prediction.score.max)
atac$predicted.id <- factor(atac$predicted.id) # to make the colors matchData integration plots in Seurat space
pat1 <- DimPlot(atac, reduction = "umap", group.by="treatment", label=TRUE) +
NoLegend() +
ggtitle("scATAC-seq truth") +
xlim(c(-6,7))
pat2 <- DimPlot(atac, group.by = "predicted.id", label = TRUE, repel = TRUE) +
ggtitle("scATAC-seq predicted") +
NoLegend() +
scale_colour_hue(drop = FALSE) +
xlim(c(-6,7))
pat3 <- DimPlot(atac, group.by = "clHBCMerged", label = TRUE, repel = TRUE) +
ggtitle("scATAC-seq ArchR clusters") +
NoLegend() +
scale_colour_hue(drop = FALSE) +
xlim(c(-6,7))
prna4 <- DimPlot(rna, group.by = "celltype", label = TRUE, repel = TRUE) + ggtitle("scRNA-seq") +
NoLegend()
pat1 + pat2 + pat3 + prna4## Warning: Removed 14 rows containing missing values (geom_point).
## Warning: Removed 14 rows containing missing values (geom_point).
## Warning: Removed 14 rows containing missing values (geom_point).
Data integration plots in ArchR space
# prediction in ArchR DR space
archrUmap$predicted.id <- celltype.predictions$predicted.id
pat1_2 <- ggplot(archrUmap, aes(x=UMAP1, y=UMAP2, color=treatment)) +
geom_point(size = .75) +
ggtitle("scATAC-seq treatment") +
theme_classic() +
scale_color_manual(values = c("orange", "darkseagreen3"))
pat2_2 <- ggplot(archrUmap, aes(x=UMAP1, y=UMAP2, color=predicted.id)) +
geom_point(size = .75) +
ggtitle("scATAC-seq predicted") +
theme_classic()
pat3_2 <- ggplot(archrUmap, aes(x=UMAP1, y=UMAP2, color=clHBCMerged)) +
geom_point(size = .75) +
ggtitle("scATAC-seq clusters") +
theme_classic()
prna1_2 <- DimPlot(rna, group.by = "celltype", label = TRUE, repel = TRUE) +
ggtitle("scRNA-seq") +
NoLegend()
pat1_2 + pat2_2 + pat3_2 + prna1_2
Joint dimensionality reduction
genes.use <- VariableFeatures(rna)
refdata <- GetAssayData(rna, assay = "RNA", slot = "data")[genes.use, ]
# refdata (input) contains a scRNA-seq expression matrix for the scRNA-seq cells. imputation
# (output) will contain an imputed scRNA-seq matrix for each of the ATAC cells
imputation <- TransferData(anchorset = transfer.anchors,
refdata = refdata,
weight.reduction = atac[["lsi"]])## Finding integration vectors## Finding integration vector weights## Transfering 2000 features onto reference data# this line adds the imputed data matrix to the pbmc.atac object
atacJoint <- atac
atacJoint[["RNA"]] <- imputation
coembed <- merge(x = rna, y = atacJoint)
# Finally, we run PCA and UMAP on this combined object, to visualize the co-embedding of both
# datasets
coembed <- ScaleData(coembed, features = genes.use, do.scale = FALSE)## Centering data matrixcoembed <- RunPCA(coembed, features = genes.use, verbose = FALSE)
coembed <- RunUMAP(coembed, dims = 1:30)## 15:13:47 UMAP embedding parameters a = 0.9922 b = 1.112## 15:13:47 Read 7140 rows and found 30 numeric columns## 15:13:47 Using Annoy for neighbor search, n_neighbors = 30## 15:13:47 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 15:13:49 Writing NN index file to temp file /var/folders/15/gs2v829x4gxcdv0tlgk_5l4m0000gn/T//RtmpPH7YMA/file13a0979f7df52
## 15:13:49 Searching Annoy index using 1 thread, search_k = 3000
## 15:13:51 Annoy recall = 100%
## 15:13:52 Commencing smooth kNN distance calibration using 1 thread
## 15:13:54 Initializing from normalized Laplacian + noise
## 15:13:54 Commencing optimization for 500 epochs, with 318608 positive edges
## 15:14:05 Optimization finishedcoembed$celltype <- ifelse(!is.na(coembed$celltype), coembed$celltype, coembed$predicted.id)
p1 <- DimPlot(coembed, group.by = "tech")
p2 <- DimPlot(coembed, group.by = "celltype", label = TRUE, repel = TRUE)
p1 + p2
Compare predictions with truth and ArchR clusters
library(ggalluvial)
dfCellWide <- data.frame(treatment=atac$treatment,
predicted=atac$predicted.id,
cluster= archProj@cellColData$clHBCMerged)
pAlluv <- ggplot(dfCellWide,
aes( axis1 = cluster, axis2 = predicted, axis3=treatment)) +
geom_alluvium(aes(fill=treatment), width = 1/12) +
scale_x_discrete(limits = c("Cluster", "Predicted", "Treatment"), expand = c(.05, .05)) +
geom_stratum(width = 1/12, fill = "white", color = "grey") +
geom_label(stat = "stratum", aes(label = after_stat(stratum))) +
theme(panel.background = element_rect(fill="white"))
pAlluv
Check gene expression of chromatin-responsive genes upon injury
Here, we check the gene expression for genes that were found to be DA in the scATAC-seq data when testing for the injury effect for both the injured vs uninjured cluster and injury vs uninjured cells in the hybrid cluster.
Chromatin more accessible in injury
chromUpInjury <- read.table("../data/scATAC/sharedUpInjury.txt",
stringsAsFactors = FALSE)[,1]
chromUpInjury <- chromUpInjury[chromUpInjury %in% rownames(countHBC)]
plotByExpr <- function(rna, genes, counts, name="feature", scale=TRUE){
genes <- genes[genes %in% rownames(countHBC)]
if(scale){
scaledCounts <- t(scale(t(counts[genes,]+1e-5)))
scaledCounts <- scaledCounts[!is.na(rowVars(scaledCounts)),]
yhat <- colSums(scaledCounts)
} else {
yhat <- colSums(counts[genes,]+1e-5)
}
rna <- AddMetaData(rna, yhat, col.name = name)
FeaturePlot(rna, features = name) + theme(legend.position = "none")
}
plotByExpr(rna, chromUpInjury, countHBC, "chromUpInjurySum", scale=FALSE)
plotByExpr(rna, chromUpInjury, countHBC, "chromUpInjurySum_scaled", scale=TRUE)
Chromatin less accessible in injury
chromDownInjury <- read.table("../data/scATAC/sharedDownInjury.txt",
stringsAsFactors = FALSE)[,1]
chromDownInjury <- chromDownInjury[chromDownInjury %in% rownames(countHBC)]
plotByExpr(rna, chromDownInjury, countHBC, "chromDownInjurySum", scale=FALSE)
plotByExpr(rna, chromDownInjury, countHBC, "chromDownjurySum_scaled", scale=TRUE)
Check gene expression of chromatin-responsive genes in hybrid vs uninjured ATAC-seq clusters
chromHybrid <- read.table("../data/scATAC/sharedHybridGenes.txt",
stringsAsFactors = FALSE)[,1]
chromHybrid <- chromHybrid[chromHybrid %in% rownames(countHBC)]
plotByExpr(rna, chromHybrid, countHBC, "hybridSum", scale=FALSE)
plotByExpr(rna, chromHybrid, countHBC, "chybridSum_scaled", scale=TRUE)
ATAC-seq cluster markers
# atacMarkers <- readRDS("../../../scATAC/data/markerList_atac_ArchR_hbcClusters.rds")
atacMarkers <- readRDS("../data/scATAC/markerList_atac_ArchR_hbcClusters.rds")
pc1 <- plotByExpr(rna, atacMarkers$C1_Inj$name, countHBC, "C1_Inj", scale=TRUE)
pc2 <- plotByExpr(rna, atacMarkers$C2_Hybrid$name, countHBC, "C2_Hybrid", scale=TRUE)
pc3 <- plotByExpr(rna, atacMarkers$C3_UI$name, countHBC, "C3_UI", scale=TRUE)
pMarkersATAC <- cowplot::plot_grid(plotlist=list(pc1,pc2,pc3), nrow=1)
# ggsave("../plots/ATACMarkersOnRNA.png", width=19)
## are the ATAC markers driven by both injured and uninjured cells in hybrid cluster? Yes.
activityHybridMarkers <- as.matrix(atac[["ACTIVITY"]]@counts)[atacMarkers$C2_Hybrid$name,]
df <- data.frame(y=colSums(t(scale(t(activityHybridMarkers)))),
state=factor(interaction(atac$clHBCMerged, atac$treatment)))
ggplot(df, aes(x=state, y=y)) +
geom_boxplot() +
xlab("Group (cell state x treatment)") +
ylab("Scaled activity of hybrid markers") +
ggtitle("Cell state expression of hybrid markers") +
theme_bw() +
theme(axis.text.x = element_text(angle=45, vjust=0.5))
# ggsave("../plots/atacHybridMarkerCellContribution.png", width=9)Redo cell label prediction when assigning hybrid labels
For reference, old plot
# Seurat
pat1 + pat2 + pat3 + prna4## Warning: Removed 14 rows containing missing values (geom_point).
## Warning: Removed 14 rows containing missing values (geom_point).
## Warning: Removed 14 rows containing missing values (geom_point).
# ArchR DR space
pat1_2 + pat2_2 + pat3_2 + prna1_2
Per hybrid cluster separately
ctHybridSub <- as.character(rna$celltype)
hlpid <- which(rna@reductions$umap@cell.embeddings[,1] > 2 & rna@reductions$umap@cell.embeddings[,2]>-1)
ctHybridSub[hlpid] <- NA
hlp2id <- which(rna@reductions$umap@cell.embeddings[,1] > 2 & rna@reductions$umap@cell.embeddings[,2]>-1 & rna@reductions$umap@cell.embeddings[,2]<4)
ctHybridSub[hlp2id] <- "hybrid1"
hlp1id <- which(rna@reductions$umap@cell.embeddings[,1] > 2 & rna@reductions$umap@cell.embeddings[,2]>4)
ctHybridSub[hlp1id] <- "hybrid2"
rna$ctHybridSub <- ctHybridSub
# ctHybridSub <- as.character(rna$celltype)
# hlpid <- which(rna@reductions$umap@cell.embeddings[,1] > 2 & rna@reductions$umap@cell.embeddings[,2] < 2)
# ctHybridSub[hlpid] <- NA
# hlp2id <- which(rna@reductions$umap@cell.embeddings[,1] > 2 & rna@reductions$umap@cell.embeddings[,2] > -1 & rna@reductions$umap@cell.embeddings[,2]<2)
# ctHybridSub[hlp2id] <- "hybrid1"
# hlp1id <- which(rna@reductions$umap@cell.embeddings[,1] > 2 & rna@reductions$umap@cell.embeddings[,2] < -4)
# ctHybridSub[hlp1id] <- "hybrid2"
# rna$ctHybridSub <- ctHybridSub
DimPlot(rna, group.by = "ctHybridSub", label = FALSE, repel = TRUE) +
ggtitle("scRNA-seq") 
# predict ATAC cell type based on RNA labels
celltype.predictions3 <- TransferData(anchorset = transfer.anchors,
refdata = rna$ctHybridSub,
weight.reduction = atac[["lsi"]])## Finding integration vectors## Finding integration vector weights## Predicting cell labels# add predicted labels to ATAC
atac$predictedCtHybrid <- celltype.predictions3$predicted.id
atac$predScore <- celltype.predictions3$prediction.score.max
atac$predicted.id <- factor(atac$predicted.id) # to make the colors match
# prediction in Seurat RD space
ppred1 <- DimPlot(atac, reduction = "umap", group.by="treatment", label=FALSE) +
ggtitle("scATAC-seq treatment") +
xlim(c(-6, 7))
ppred2 <- DimPlot(atac, group.by = "predictedCtHybrid", label = FALSE, repel = TRUE) +
ggtitle("scATAC-seq predicted state") +
scale_colour_hue(drop = FALSE) +
xlim(c(-6, 7))
ppred3 <- DimPlot(atac, group.by = "clHBCMerged", label = FALSE, repel = TRUE) +
ggtitle("scATAC-seq cell state") +
scale_colour_hue(drop = FALSE) +
xlim(c(-6, 7))
ppred4 <- DimPlot(rna, group.by = "ctHybridSub", label = FALSE, repel = TRUE) +
ggtitle("scRNA-seq")
ppred1 + ppred2 + ppred3 + ppred4## Warning: Removed 14 rows containing missing values (geom_point).
## Warning: Removed 14 rows containing missing values (geom_point).
## Warning: Removed 14 rows containing missing values (geom_point).
# ggsave("../plots/cellTransferSubHybrid.png", width=10)
## size of cell is prediction prob
DimPlot(atac, group.by = "predictedCtHybrid", label = TRUE, repel = TRUE, pt.size=celltype.predictions3$prediction.score.max*1.25) + ggtitle("scATAC-seq predicted") +
NoLegend() + scale_colour_hue(drop = FALSE) + xlim(c(-6,6))## Warning: Removed 55 rows containing missing values (geom_point).
boxplot(celltype.predictions3$prediction.score.max ~ celltype.predictions3$predicted.id,
xlab="Predicted state", ylab="Probability")
# prediction in ArchR DR space
archrUmap$predicted.id <- celltype.predictions3$predicted.id
ppred1_2 <- ggplot(archrUmap, aes(x=UMAP1, y=UMAP2, color=treatment)) +
geom_point(size = .75) +
ggtitle("scATAC-seq treatment") +
theme_classic() +
scale_color_manual(values = c("orange", "darkseagreen3"))
ppred2_2 <- ggplot(archrUmap, aes(x=UMAP1, y=UMAP2, color=predicted.id)) +
geom_point(size = .5) +
ggtitle("scATAC-seq predicted") +
theme_classic()
ppred3_2 <- ggplot(archrUmap, aes(x=UMAP1, y=UMAP2, color=clHBCMerged)) +
geom_point(size = .75) +
ggtitle("scATAC-seq clusters") +
theme_classic()
ppred4_2 <- DimPlot(rna, group.by = "ctHybridSub", label = FALSE, repel = TRUE, pt.size=.1) +
ggtitle("scRNA-seq") +
theme(axis.title = element_text(size = 11),
axis.text = element_text(size = 9),
plot.title = element_text(face = "plain", size=13))
ppred1_2 + ppred2_2 + ppred3_2 + ppred4_2
ggsave("../plots/predictionsRNAATAC_archRSpace.png", width=12, height=9)
## size of cell is prediction prob
ggplot(archrUmap, aes(x=UMAP1, y=UMAP2, color=predicted.id)) +
geom_point(size = celltype.predictions3$prediction.score.max*1.25) +
ggtitle("scATAC-seq predicted") +
theme_classic()
boxplot(celltype.predictions3$prediction.score.max ~ celltype.predictions3$predicted.id,
xlab="Predicted state", ylab="Probability")
## distribution of injured and uninjured cells in each cluster
df <- data.frame(treatment=atac$treatment,
cluster=atac$predictedCtHybrid)
dfTotal <- df %>%
group_by(cluster) %>%
mutate(totalCells=n()) %>%
group_by(cluster, treatment) %>%
summarize(fracCells = n() / unique(totalCells))## `summarise()` has grouped output by 'cluster'. You can override using the `.groups` argument.ggplot(dfTotal, aes(x=treatment, fill=treatment, y=fracCells)) +
geom_bar(position="dodge", stat="identity") +
scale_fill_manual(values = c("orange", "darkseagreen3")) +
facet_wrap(.~cluster) +
ylab("Fraction of cells") +
theme_classic()
Joint dimensionality reduction
atacJoint <- atac
DefaultAssay(atacJoint) <- "ATAC"
VariableFeatures(atacJoint) <- names(which(Matrix::rowSums(atacJoint) > 100))
atacJoint <- RunLSI(atacJoint, n = 50, scale.max = NULL)## Warning: RunLSI is being moved to Signac. Equivalent functionality can be
## achieved via the Signac::RunTFIDF and Signac::RunSVD functions; for more
## information on Signac, please see https://github.com/timoast/Signac
## Warning: RunLSI is being moved to Signac. Equivalent functionality can be
## achieved via the Signac::RunTFIDF and Signac::RunSVD functions; for more
## information on Signac, please see https://github.com/timoast/Signac
## Warning: RunLSI is being moved to Signac. Equivalent functionality can be
## achieved via the Signac::RunTFIDF and Signac::RunSVD functions; for more
## information on Signac, please see https://github.com/timoast/Signac## Performing TF-IDF normalization## Running SVD on TF-IDF matrix## Scaling cell embeddingsatacJoint <- RunUMAP(atacJoint, reduction = "lsi", dims = 2:50)## 15:18:08 UMAP embedding parameters a = 0.9922 b = 1.112## 15:18:08 Read 4076 rows and found 49 numeric columns## 15:18:08 Using Annoy for neighbor search, n_neighbors = 30## 15:18:08 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 15:18:09 Writing NN index file to temp file /var/folders/15/gs2v829x4gxcdv0tlgk_5l4m0000gn/T//RtmpPH7YMA/file13a09478fc41c
## 15:18:09 Searching Annoy index using 1 thread, search_k = 3000
## 15:18:10 Annoy recall = 100%
## 15:18:13 Commencing smooth kNN distance calibration using 1 thread
## 15:18:15 Initializing from normalized Laplacian + noise
## 15:18:15 Commencing optimization for 500 epochs, with 155968 positive edges
## 15:18:21 Optimization finishedgenes.use <- VariableFeatures(rna)
refdata <- GetAssayData(rna, assay = "RNA", slot = "data")[genes.use, ]
# refdata (input) contains a scRNA-seq expression matrix for the scRNA-seq cells. imputation
# (output) will contain an imputed scRNA-seq matrix for each of the ATAC cells
imputation <- TransferData(anchorset = transfer.anchors,
refdata = refdata,
weight.reduction = atacJoint[["lsi"]])## Finding integration vectors
## Finding integration vector weights
## Transfering 2000 features onto reference data# this line adds the imputed data matrix to the pbmc.atac object
atacJoint[["RNA"]] <- imputation
coembed <- merge(x = rna, y = atacJoint)
# Finally, we run PCA and UMAP on this combined object, to visualize the co-embedding of both
# datasets
coembed <- ScaleData(coembed, features = genes.use, do.scale = FALSE)## Centering data matrixcoembed <- RunPCA(coembed, features = genes.use, verbose = FALSE)
coembed <- RunUMAP(coembed, dims = 1:30)## 15:18:58 UMAP embedding parameters a = 0.9922 b = 1.112
## 15:18:58 Read 7140 rows and found 30 numeric columns
## 15:18:58 Using Annoy for neighbor search, n_neighbors = 30
## 15:18:58 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 15:18:59 Writing NN index file to temp file /var/folders/15/gs2v829x4gxcdv0tlgk_5l4m0000gn/T//RtmpPH7YMA/file13a092f942ef8
## 15:18:59 Searching Annoy index using 1 thread, search_k = 3000
## 15:19:01 Annoy recall = 100%
## 15:19:03 Commencing smooth kNN distance calibration using 1 thread
## 15:19:04 Initializing from normalized Laplacian + noise
## 15:19:04 Commencing optimization for 500 epochs, with 323034 positive edges
## 15:19:16 Optimization finishedcelltype.predictionsJoint <- TransferData(anchorset = transfer.anchors,
refdata = rna$ctHybridSub,
weight.reduction = atacJoint[["lsi"]],
dims=1:50)## Finding integration vectors
## Finding integration vector weights
## Predicting cell labels# fill in predicted labels: using joint DR or procedure above is very similar.
coembed$ctHybridSubPred <- coembed$ctHybridSub
## based on joint DR
# coembed@meta.data[rownames(celltype.predictionsJoint), "ctHybridSubPred"] <- celltype.predictionsJoint$predicted.id
## ones we had before
coembed@meta.data[rownames(celltype.predictions3), "ctHybridSubPred"] <- celltype.predictions3$predicted.id
p1 <- DimPlot(coembed, group.by = "tech") + ggtitle("Modality")
p2 <- DimPlot(coembed, group.by = "ctHybridSub", label = TRUE, repel = TRUE) + ggtitle("RNA clusters")
p3 <- DimPlot(coembed, group.by = "clHBCMerged", label = TRUE, repel = TRUE) + ggtitle("ATAC clusters")
p4 <- DimPlot(coembed, group.by = "ctHybridSubPred", label = TRUE, repel = TRUE) + ggtitle("Transferred clusters")
p1 + p2 + p3 + p4
ggsave("../plots/jointDimRed_v1.png", width=10, height=6)Joint dimensionality reduction using LIGER
# liger basically takes two count matrices as input.
# we'll use the rounded ACTIVITY count matrix as input for scATAC.
atacCounts <- round(atac[["ACTIVITY"]]@counts)
rnaCounts <- round(rna@assays$RNA@counts)
library(liger)##
## Attaching package: 'liger'## The following objects are masked from 'package:scater':
##
## normalize, runTSNE, runUMAP## The following object is masked from 'package:BiocGenerics':
##
## normalize# hbcData <- list(atac = atacCounts[rownames(atacCounts) %in% rownames(rnaCounts),],
# rna = rnaCounts)
# hbcData <- createLiger(hbcData)
# # preprocess
# hbcData <- normalize(hbcData)
# hbcData <- selectGenes(hbcData, datasets.use = 2)
# hbcData <- scaleNotCenter(hbcData)
# scaleNotCenter takes extremely long. Result is therefore saved.
hbcData <- readRDS("../data/hbcData_liger_afterScaleNotCenter.rds")
# suggestLambda(hbcData, k=30)
hbcData <- optimizeALS(hbcData, k=30, lambda = 10)##
|
| | 0%
|
|== | 3%
|
|===== | 7%
|
|======= | 10%
|
|========= | 13%
|
|============ | 17%
|
|============== | 20%
|
|================ | 23%
|
|=================== | 27%
|
|===================== | 30%
|
|======================= | 33%
|
|========================== | 37%
|
|============================ | 40%
|
|============================== | 43%
|
|================================= | 47%
|
|=================================== | 50%
|
|===================================== | 53%
|
|======================================== | 57%
|
|========================================== | 60%
|
|============================================ | 63%
|
|=============================================== | 67%
|
|================================================= | 70%
|
|=================================================== | 73%
|
|====================================================== | 77%
|
|======================================================== | 80%
|
|========================================================== | 83%
|
|============================================================= | 87%
|
|=============================================================== | 90%
|
|================================================================= | 93%
|
|==================================================================== | 97%
|
|======================================================================| 100%
## Finished in 9.718121 mins, 30 iterations.
## Max iterations set: 30.
## Final objective delta: 3.445833e-05.
## Best results with seed 1.hbcData <- quantile_norm(hbcData)
hbcData <- runUMAP(hbcData, distance='cosine', min_dist=.3)
clSeurat <- coembed@meta.data$ctHybridSubPred
names(clSeurat) <- rownames(coembed@meta.data)
clSeurat <- clSeurat[names(hbcData@clusters)]
liger::plotByDatasetAndCluster(hbcData, clusters = clSeurat)
Markers for sub-hybrid clusters
## add annotation according to prediction
predictedAnn <- archProj@cellColData$clHBCMerged
predictedAnn[predictedAnn == "C3_UI"] <- "resting"
predictedAnn[predictedAnn == "C1_inj"] <- "activated"
hybId <- predictedAnn == "C2_hybrid"
predictedAnn[which(hybId)] <- NA
predictedAnn[hybId &
!celltype.predictions3$predicted.id == "hybrid1" &
!celltype.predictions3$predicted.id == "hybrid2"] <- "hybrid"
predictedAnn[which(celltype.predictions3$predicted.id == "hybrid1")] <- "hybrid1"
predictedAnn[which(celltype.predictions3$predicted.id == "hybrid2")] <- "hybrid2"
table(predictedAnn)## predictedAnn
## C1_Inj C2_Hybrid hybrid1 hybrid2 resting
## 1485 933 197 223 1238## visualize with Seurat
atac$predAnn <- predictedAnn
DimPlot(atac, group.by = "predAnn", label = TRUE, repel = TRUE) + ggtitle("scATAC-seq predicted") +
NoLegend() + scale_colour_hue(drop = FALSE) + xlim(c(-6,7))## Warning: Removed 14 rows containing missing values (geom_point).
## distribution of injured and uninjured cells in each cluster
df <- data.frame(treatment=atac$treatment,
cluster=atac$predAnn)
dfTotal <- df %>%
group_by(cluster) %>%
mutate(totalCells=n()) %>%
group_by(cluster, treatment) %>%
summarize(fracCells = n() / unique(totalCells))## `summarise()` has grouped output by 'cluster'. You can override using the `.groups` argument.ggplot(dfTotal, aes(x=treatment, fill=treatment, y=fracCells)) +
geom_bar(position="dodge", stat="identity") +
facet_wrap(.~cluster) +
ylab("Fraction of cells")
# saveRDS(celltype.predictions3, file="../data/scATAC/20201105_celltypePredictions3.rds")
## are the ATAC markers driven by both injured and uninjured cells in all hybrid cluster? Yes.
activityHybridMarkers <- as.matrix(atac[["ACTIVITY"]]@counts)[atacMarkers$C2_Hybrid$name,]
df <- data.frame(y=colSums(t(scale(t(activityHybridMarkers)))),
state=factor(interaction(atac$predAnn, atac$treatment)))
ggplot(df, aes(x=state, y=y)) +
geom_boxplot() +
xlab("Group (cell state x treatment)") +
ylab("Scaled activity of hybrid markers") +
ggtitle("Cell state expression of hybrid markers") +
theme_bw() +
theme(axis.text.x = element_text(angle=45, vjust=0.5))
# ggsave("../plots/atacHybridMarkerCellContributionAllHybridClusters.png", width=9)New peak calling based on sub-clustering of hybrid cells as suggested by RNA
archProj@cellColData$predAnn <- predictedAnn
# pseudobulk for peak calling
archProj <- addGroupCoverages(ArchRProj = archProj, groupBy = "predAnn",
minReplicates = 2, maxReplicates = 6,
force = TRUE,
verbose = FALSE)## ArchR logging to : ArchRLogs/ArchR-addGroupCoverages-13a09da455d4-Date-2021-04-02_Time-15-29-29.log
## If there is an issue, please report to github with logFile!## C1_Inj (1 of 5) : CellGroups N = 4## C2_Hybrid (2 of 5) : CellGroups N = 6## hybrid1 (3 of 5) : CellGroups N = 3## hybrid2 (4 of 5) : CellGroups N = 2## resting (5 of 5) : CellGroups N = 5## 2021-04-02 15:29:30 : Creating Coverage Files!, 0.022 mins elapsed.## 2021-04-02 15:29:30 : Batch Execution w/ safelapply!, 0.022 mins elapsed.## Number of Cells = 500## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 263## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 163## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 54## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 269## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 252## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 140## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 113## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 83## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 76## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 94## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 57## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 46## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 170## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 53## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 500## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 355## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 178## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 54## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 51## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## 2021-04-02 15:36:46 : Adding Kmer Bias to Coverage Files!, 7.28 mins elapsed.## Kmer Bias chr1 (1 of 21)## chr1 Coverage File chr1 (1 of 20)
## Coverage File chr1 (2 of 20)
## Coverage File chr1 (3 of 20)
## Coverage File chr1 (4 of 20)
## Coverage File chr1 (5 of 20)
## Coverage File chr1 (6 of 20)
## Coverage File chr1 (7 of 20)
## Coverage File chr1 (8 of 20)
## Coverage File chr1 (9 of 20)
## Coverage File chr1 (10 of 20)
## Coverage File chr1 (11 of 20)
## Coverage File chr1 (12 of 20)
## Coverage File chr1 (13 of 20)
## Coverage File chr1 (14 of 20)
## Coverage File chr1 (15 of 20)
## Coverage File chr1 (16 of 20)
## Coverage File chr1 (17 of 20)
## Coverage File chr1 (18 of 20)
## Coverage File chr1 (19 of 20)
## Coverage File chr1 (20 of 20)
## Kmer Bias chr10 (2 of 21)
## chr10 Coverage File chr10 (1 of 20)
## Coverage File chr10 (2 of 20)
## Coverage File chr10 (3 of 20)
## Coverage File chr10 (4 of 20)
## Coverage File chr10 (5 of 20)
## Coverage File chr10 (6 of 20)
## Coverage File chr10 (7 of 20)
## Coverage File chr10 (8 of 20)
## Coverage File chr10 (9 of 20)
## Coverage File chr10 (10 of 20)
## Coverage File chr10 (11 of 20)
## Coverage File chr10 (12 of 20)
## Coverage File chr10 (13 of 20)
## Coverage File chr10 (14 of 20)
## Coverage File chr10 (15 of 20)
## Coverage File chr10 (16 of 20)
## Coverage File chr10 (17 of 20)
## Coverage File chr10 (18 of 20)
## Coverage File chr10 (19 of 20)
## Coverage File chr10 (20 of 20)
## Kmer Bias chr11 (3 of 21)
## chr11 Coverage File chr11 (1 of 20)
## Coverage File chr11 (2 of 20)
## Coverage File chr11 (3 of 20)
## Coverage File chr11 (4 of 20)
## Coverage File chr11 (5 of 20)
## Coverage File chr11 (6 of 20)
## Coverage File chr11 (7 of 20)
## Coverage File chr11 (8 of 20)
## Coverage File chr11 (9 of 20)
## Coverage File chr11 (10 of 20)
## Coverage File chr11 (11 of 20)
## Coverage File chr11 (12 of 20)
## Coverage File chr11 (13 of 20)
## Coverage File chr11 (14 of 20)
## Coverage File chr11 (15 of 20)
## Coverage File chr11 (16 of 20)
## Coverage File chr11 (17 of 20)
## Coverage File chr11 (18 of 20)
## Coverage File chr11 (19 of 20)
## Coverage File chr11 (20 of 20)
## Kmer Bias chr12 (4 of 21)
## chr12 Coverage File chr12 (1 of 20)
## Coverage File chr12 (2 of 20)
## Coverage File chr12 (3 of 20)
## Coverage File chr12 (4 of 20)
## Coverage File chr12 (5 of 20)
## Coverage File chr12 (6 of 20)
## Coverage File chr12 (7 of 20)
## Coverage File chr12 (8 of 20)
## Coverage File chr12 (9 of 20)
## Coverage File chr12 (10 of 20)
## Coverage File chr12 (11 of 20)
## Coverage File chr12 (12 of 20)
## Coverage File chr12 (13 of 20)
## Coverage File chr12 (14 of 20)
## Coverage File chr12 (15 of 20)
## Coverage File chr12 (16 of 20)
## Coverage File chr12 (17 of 20)
## Coverage File chr12 (18 of 20)
## Coverage File chr12 (19 of 20)
## Coverage File chr12 (20 of 20)
## Kmer Bias chr13 (5 of 21)
## chr13 Coverage File chr13 (1 of 20)
## Coverage File chr13 (2 of 20)
## Coverage File chr13 (3 of 20)
## Coverage File chr13 (4 of 20)
## Coverage File chr13 (5 of 20)
## Coverage File chr13 (6 of 20)
## Coverage File chr13 (7 of 20)
## Coverage File chr13 (8 of 20)
## Coverage File chr13 (9 of 20)
## Coverage File chr13 (10 of 20)
## Coverage File chr13 (11 of 20)
## Coverage File chr13 (12 of 20)
## Coverage File chr13 (13 of 20)
## Coverage File chr13 (14 of 20)
## Coverage File chr13 (15 of 20)
## Coverage File chr13 (16 of 20)
## Coverage File chr13 (17 of 20)
## Coverage File chr13 (18 of 20)
## Coverage File chr13 (19 of 20)
## Coverage File chr13 (20 of 20)
## Kmer Bias chr14 (6 of 21)
## chr14 Coverage File chr14 (1 of 20)
## Coverage File chr14 (2 of 20)
## Coverage File chr14 (3 of 20)
## Coverage File chr14 (4 of 20)
## Coverage File chr14 (5 of 20)
## Coverage File chr14 (6 of 20)
## Coverage File chr14 (7 of 20)
## Coverage File chr14 (8 of 20)
## Coverage File chr14 (9 of 20)
## Coverage File chr14 (10 of 20)
## Coverage File chr14 (11 of 20)
## Coverage File chr14 (12 of 20)
## Coverage File chr14 (13 of 20)
## Coverage File chr14 (14 of 20)
## Coverage File chr14 (15 of 20)
## Coverage File chr14 (16 of 20)
## Coverage File chr14 (17 of 20)
## Coverage File chr14 (18 of 20)
## Coverage File chr14 (19 of 20)
## Coverage File chr14 (20 of 20)
## Kmer Bias chr15 (7 of 21)
## chr15 Coverage File chr15 (1 of 20)
## Coverage File chr15 (2 of 20)
## Coverage File chr15 (3 of 20)
## Coverage File chr15 (4 of 20)
## Coverage File chr15 (5 of 20)
## Coverage File chr15 (6 of 20)
## Coverage File chr15 (7 of 20)
## Coverage File chr15 (8 of 20)
## Coverage File chr15 (9 of 20)
## Coverage File chr15 (10 of 20)
## Coverage File chr15 (11 of 20)
## Coverage File chr15 (12 of 20)
## Coverage File chr15 (13 of 20)
## Coverage File chr15 (14 of 20)
## Coverage File chr15 (15 of 20)
## Coverage File chr15 (16 of 20)
## Coverage File chr15 (17 of 20)
## Coverage File chr15 (18 of 20)
## Coverage File chr15 (19 of 20)
## Coverage File chr15 (20 of 20)
## Kmer Bias chr16 (8 of 21)
## chr16 Coverage File chr16 (1 of 20)
## Coverage File chr16 (2 of 20)
## Coverage File chr16 (3 of 20)
## Coverage File chr16 (4 of 20)
## Coverage File chr16 (5 of 20)
## Coverage File chr16 (6 of 20)
## Coverage File chr16 (7 of 20)
## Coverage File chr16 (8 of 20)
## Coverage File chr16 (9 of 20)
## Coverage File chr16 (10 of 20)
## Coverage File chr16 (11 of 20)
## Coverage File chr16 (12 of 20)
## Coverage File chr16 (13 of 20)
## Coverage File chr16 (14 of 20)
## Coverage File chr16 (15 of 20)
## Coverage File chr16 (16 of 20)
## Coverage File chr16 (17 of 20)
## Coverage File chr16 (18 of 20)
## Coverage File chr16 (19 of 20)
## Coverage File chr16 (20 of 20)
## Kmer Bias chr17 (9 of 21)
## chr17 Coverage File chr17 (1 of 20)
## Coverage File chr17 (2 of 20)
## Coverage File chr17 (3 of 20)
## Coverage File chr17 (4 of 20)
## Coverage File chr17 (5 of 20)
## Coverage File chr17 (6 of 20)
## Coverage File chr17 (7 of 20)
## Coverage File chr17 (8 of 20)
## Coverage File chr17 (9 of 20)
## Coverage File chr17 (10 of 20)
## Coverage File chr17 (11 of 20)
## Coverage File chr17 (12 of 20)
## Coverage File chr17 (13 of 20)
## Coverage File chr17 (14 of 20)
## Coverage File chr17 (15 of 20)
## Coverage File chr17 (16 of 20)
## Coverage File chr17 (17 of 20)
## Coverage File chr17 (18 of 20)
## Coverage File chr17 (19 of 20)
## Coverage File chr17 (20 of 20)
## Kmer Bias chr18 (10 of 21)
## chr18 Coverage File chr18 (1 of 20)
## Coverage File chr18 (2 of 20)
## Coverage File chr18 (3 of 20)
## Coverage File chr18 (4 of 20)
## Coverage File chr18 (5 of 20)
## Coverage File chr18 (6 of 20)
## Coverage File chr18 (7 of 20)
## Coverage File chr18 (8 of 20)
## Coverage File chr18 (9 of 20)
## Coverage File chr18 (10 of 20)
## Coverage File chr18 (11 of 20)
## Coverage File chr18 (12 of 20)
## Coverage File chr18 (13 of 20)
## Coverage File chr18 (14 of 20)
## Coverage File chr18 (15 of 20)
## Coverage File chr18 (16 of 20)
## Coverage File chr18 (17 of 20)
## Coverage File chr18 (18 of 20)
## Coverage File chr18 (19 of 20)
## Coverage File chr18 (20 of 20)
## Kmer Bias chr19 (11 of 21)
## chr19 Coverage File chr19 (1 of 20)
## Coverage File chr19 (2 of 20)
## Coverage File chr19 (3 of 20)
## Coverage File chr19 (4 of 20)
## Coverage File chr19 (5 of 20)
## Coverage File chr19 (6 of 20)
## Coverage File chr19 (7 of 20)
## Coverage File chr19 (8 of 20)
## Coverage File chr19 (9 of 20)
## Coverage File chr19 (10 of 20)
## Coverage File chr19 (11 of 20)
## Coverage File chr19 (12 of 20)
## Coverage File chr19 (13 of 20)
## Coverage File chr19 (14 of 20)
## Coverage File chr19 (15 of 20)
## Coverage File chr19 (16 of 20)
## Coverage File chr19 (17 of 20)
## Coverage File chr19 (18 of 20)
## Coverage File chr19 (19 of 20)
## Coverage File chr19 (20 of 20)
## Kmer Bias chr2 (12 of 21)
## chr2 Coverage File chr2 (1 of 20)
## Coverage File chr2 (2 of 20)
## Coverage File chr2 (3 of 20)
## Coverage File chr2 (4 of 20)
## Coverage File chr2 (5 of 20)
## Coverage File chr2 (6 of 20)
## Coverage File chr2 (7 of 20)
## Coverage File chr2 (8 of 20)
## Coverage File chr2 (9 of 20)
## Coverage File chr2 (10 of 20)
## Coverage File chr2 (11 of 20)
## Coverage File chr2 (12 of 20)
## Coverage File chr2 (13 of 20)
## Coverage File chr2 (14 of 20)
## Coverage File chr2 (15 of 20)
## Coverage File chr2 (16 of 20)
## Coverage File chr2 (17 of 20)
## Coverage File chr2 (18 of 20)
## Coverage File chr2 (19 of 20)
## Coverage File chr2 (20 of 20)
## Kmer Bias chr3 (13 of 21)
## chr3 Coverage File chr3 (1 of 20)
## Coverage File chr3 (2 of 20)
## Coverage File chr3 (3 of 20)
## Coverage File chr3 (4 of 20)
## Coverage File chr3 (5 of 20)
## Coverage File chr3 (6 of 20)
## Coverage File chr3 (7 of 20)
## Coverage File chr3 (8 of 20)
## Coverage File chr3 (9 of 20)
## Coverage File chr3 (10 of 20)
## Coverage File chr3 (11 of 20)
## Coverage File chr3 (12 of 20)
## Coverage File chr3 (13 of 20)
## Coverage File chr3 (14 of 20)
## Coverage File chr3 (15 of 20)
## Coverage File chr3 (16 of 20)
## Coverage File chr3 (17 of 20)
## Coverage File chr3 (18 of 20)
## Coverage File chr3 (19 of 20)
## Coverage File chr3 (20 of 20)
## Kmer Bias chr4 (14 of 21)
## chr4 Coverage File chr4 (1 of 20)
## Coverage File chr4 (2 of 20)
## Coverage File chr4 (3 of 20)
## Coverage File chr4 (4 of 20)
## Coverage File chr4 (5 of 20)
## Coverage File chr4 (6 of 20)
## Coverage File chr4 (7 of 20)
## Coverage File chr4 (8 of 20)
## Coverage File chr4 (9 of 20)
## Coverage File chr4 (10 of 20)
## Coverage File chr4 (11 of 20)
## Coverage File chr4 (12 of 20)
## Coverage File chr4 (13 of 20)
## Coverage File chr4 (14 of 20)
## Coverage File chr4 (15 of 20)
## Coverage File chr4 (16 of 20)
## Coverage File chr4 (17 of 20)
## Coverage File chr4 (18 of 20)
## Coverage File chr4 (19 of 20)
## Coverage File chr4 (20 of 20)
## Kmer Bias chr5 (15 of 21)
## chr5 Coverage File chr5 (1 of 20)
## Coverage File chr5 (2 of 20)
## Coverage File chr5 (3 of 20)
## Coverage File chr5 (4 of 20)
## Coverage File chr5 (5 of 20)
## Coverage File chr5 (6 of 20)
## Coverage File chr5 (7 of 20)
## Coverage File chr5 (8 of 20)
## Coverage File chr5 (9 of 20)
## Coverage File chr5 (10 of 20)
## Coverage File chr5 (11 of 20)
## Coverage File chr5 (12 of 20)
## Coverage File chr5 (13 of 20)
## Coverage File chr5 (14 of 20)
## Coverage File chr5 (15 of 20)
## Coverage File chr5 (16 of 20)
## Coverage File chr5 (17 of 20)
## Coverage File chr5 (18 of 20)
## Coverage File chr5 (19 of 20)
## Coverage File chr5 (20 of 20)
## Kmer Bias chr6 (16 of 21)
## chr6 Coverage File chr6 (1 of 20)
## Coverage File chr6 (2 of 20)
## Coverage File chr6 (3 of 20)
## Coverage File chr6 (4 of 20)
## Coverage File chr6 (5 of 20)
## Coverage File chr6 (6 of 20)
## Coverage File chr6 (7 of 20)
## Coverage File chr6 (8 of 20)
## Coverage File chr6 (9 of 20)
## Coverage File chr6 (10 of 20)
## Coverage File chr6 (11 of 20)
## Coverage File chr6 (12 of 20)
## Coverage File chr6 (13 of 20)
## Coverage File chr6 (14 of 20)
## Coverage File chr6 (15 of 20)
## Coverage File chr6 (16 of 20)
## Coverage File chr6 (17 of 20)
## Coverage File chr6 (18 of 20)
## Coverage File chr6 (19 of 20)
## Coverage File chr6 (20 of 20)
## Kmer Bias chr7 (17 of 21)
## chr7 Coverage File chr7 (1 of 20)
## Coverage File chr7 (2 of 20)
## Coverage File chr7 (3 of 20)
## Coverage File chr7 (4 of 20)
## Coverage File chr7 (5 of 20)
## Coverage File chr7 (6 of 20)
## Coverage File chr7 (7 of 20)
## Coverage File chr7 (8 of 20)
## Coverage File chr7 (9 of 20)
## Coverage File chr7 (10 of 20)
## Coverage File chr7 (11 of 20)
## Coverage File chr7 (12 of 20)
## Coverage File chr7 (13 of 20)
## Coverage File chr7 (14 of 20)
## Coverage File chr7 (15 of 20)
## Coverage File chr7 (16 of 20)
## Coverage File chr7 (17 of 20)
## Coverage File chr7 (18 of 20)
## Coverage File chr7 (19 of 20)
## Coverage File chr7 (20 of 20)
## Kmer Bias chr8 (18 of 21)
## chr8 Coverage File chr8 (1 of 20)
## Coverage File chr8 (2 of 20)
## Coverage File chr8 (3 of 20)
## Coverage File chr8 (4 of 20)
## Coverage File chr8 (5 of 20)
## Coverage File chr8 (6 of 20)
## Coverage File chr8 (7 of 20)
## Coverage File chr8 (8 of 20)
## Coverage File chr8 (9 of 20)
## Coverage File chr8 (10 of 20)
## Coverage File chr8 (11 of 20)
## Coverage File chr8 (12 of 20)
## Coverage File chr8 (13 of 20)
## Coverage File chr8 (14 of 20)
## Coverage File chr8 (15 of 20)
## Coverage File chr8 (16 of 20)
## Coverage File chr8 (17 of 20)
## Coverage File chr8 (18 of 20)
## Coverage File chr8 (19 of 20)
## Coverage File chr8 (20 of 20)
## Kmer Bias chr9 (19 of 21)
## chr9 Coverage File chr9 (1 of 20)
## Coverage File chr9 (2 of 20)
## Coverage File chr9 (3 of 20)
## Coverage File chr9 (4 of 20)
## Coverage File chr9 (5 of 20)
## Coverage File chr9 (6 of 20)
## Coverage File chr9 (7 of 20)
## Coverage File chr9 (8 of 20)
## Coverage File chr9 (9 of 20)
## Coverage File chr9 (10 of 20)
## Coverage File chr9 (11 of 20)
## Coverage File chr9 (12 of 20)
## Coverage File chr9 (13 of 20)
## Coverage File chr9 (14 of 20)
## Coverage File chr9 (15 of 20)
## Coverage File chr9 (16 of 20)
## Coverage File chr9 (17 of 20)
## Coverage File chr9 (18 of 20)
## Coverage File chr9 (19 of 20)
## Coverage File chr9 (20 of 20)
## Kmer Bias chrX (20 of 21)
## chrX Coverage File chrX (1 of 20)
## Coverage File chrX (2 of 20)
## Coverage File chrX (3 of 20)
## Coverage File chrX (4 of 20)
## Coverage File chrX (5 of 20)
## Coverage File chrX (6 of 20)
## Coverage File chrX (7 of 20)
## Coverage File chrX (8 of 20)
## Coverage File chrX (9 of 20)
## Coverage File chrX (10 of 20)
## Coverage File chrX (11 of 20)
## Coverage File chrX (12 of 20)
## Coverage File chrX (13 of 20)
## Coverage File chrX (14 of 20)
## Coverage File chrX (15 of 20)
## Coverage File chrX (16 of 20)
## Coverage File chrX (17 of 20)
## Coverage File chrX (18 of 20)
## Coverage File chrX (19 of 20)
## Coverage File chrX (20 of 20)
## Kmer Bias chrY (21 of 21)
## chrY Coverage File chrY (1 of 20)
## Coverage File chrY (2 of 20)
## Coverage File chrY (3 of 20)
## Coverage File chrY (4 of 20)
## Coverage File chrY (5 of 20)
## Coverage File chrY (6 of 20)
## Coverage File chrY (7 of 20)
## Coverage File chrY (8 of 20)
## Coverage File chrY (9 of 20)
## Coverage File chrY (10 of 20)
## Coverage File chrY (11 of 20)
## Coverage File chrY (12 of 20)
## Coverage File chrY (13 of 20)
## Coverage File chrY (14 of 20)
## Coverage File chrY (15 of 20)
## Coverage File chrY (16 of 20)
## Coverage File chrY (17 of 20)
## Coverage File chrY (18 of 20)
## Coverage File chrY (19 of 20)
## Coverage File chrY (20 of 20)
## Completed Kmer Bias Calculation
## Adding Kmer Bias (1 of 20)
## Adding Kmer Bias (2 of 20)
## Adding Kmer Bias (3 of 20)
## Adding Kmer Bias (4 of 20)
## Adding Kmer Bias (5 of 20)
## Adding Kmer Bias (6 of 20)
## Adding Kmer Bias (7 of 20)
## Adding Kmer Bias (8 of 20)
## Adding Kmer Bias (9 of 20)
## Adding Kmer Bias (10 of 20)
## Adding Kmer Bias (11 of 20)
## Adding Kmer Bias (12 of 20)
## Adding Kmer Bias (13 of 20)
## Adding Kmer Bias (14 of 20)
## Adding Kmer Bias (15 of 20)
## Adding Kmer Bias (16 of 20)
## Adding Kmer Bias (17 of 20)
## Adding Kmer Bias (18 of 20)
## Adding Kmer Bias (19 of 20)
## Adding Kmer Bias (20 of 20)
## 2021-04-02 15:44:35 : Finished Creation of Coverage Files!, 15.096 mins elapsed.
## ArchR logging successful to : ArchRLogs/ArchR-addGroupCoverages-13a09da455d4-Date-2021-04-02_Time-15-29-29.log# peak calling
archProj <- addReproduciblePeakSet(
ArchRProj = archProj,
groupBy = "predAnn",
pathToMacs2 = "/Users/koenvandenberge/opt/anaconda3/bin/macs2",
verbose = FALSE)## ArchR logging to : ArchRLogs/ArchR-addReproduciblePeakSet-13a091ffed84b-Date-2021-04-02_Time-15-44-35.log
## If there is an issue, please report to github with logFile!
## Calling Peaks with Macs2
## 2021-04-02 15:44:36 : Batch Execution w/ safelapply!, 0 mins elapsed.
## 2021-04-02 15:44:36 : Group 1 of 20, Calling Peaks with MACS2!, 0 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C1_Inj._.s1Inj-1 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C1_Inj._.s1Inj-1.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-04-02 15:46:21 : Group 2 of 20, Calling Peaks with MACS2!, 1.76 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C1_Inj._.s2Inj-2 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C1_Inj._.s2Inj-2.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-04-02 15:47:53 : Group 3 of 20, Calling Peaks with MACS2!, 3.283 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C1_Inj._.s3Inj-3 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C1_Inj._.s3Inj-3.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-04-02 15:49:18 : Group 4 of 20, Calling Peaks with MACS2!, 4.712 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C1_Inj._.Other-4 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C1_Inj._.Other-4.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-04-02 15:50:12 : Group 5 of 20, Calling Peaks with MACS2!, 5.598 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C2_Hybrid._.s1Inj-5 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C2_Hybrid._.s1Inj-5.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-04-02 15:51:28 : Group 6 of 20, Calling Peaks with MACS2!, 6.876 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C2_Hybrid._.s3Inj-6 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C2_Hybrid._.s3Inj-6.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-04-02 15:52:59 : Group 7 of 20, Calling Peaks with MACS2!, 8.384 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C2_Hybrid._.s2UI-7 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C2_Hybrid._.s2UI-7.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-04-02 15:53:51 : Group 8 of 20, Calling Peaks with MACS2!, 9.262 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C2_Hybrid._.s3UI-8 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C2_Hybrid._.s3UI-8.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-04-02 15:54:39 : Group 9 of 20, Calling Peaks with MACS2!, 10.053 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C2_Hybrid._.s2Inj-9 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C2_Hybrid._.s2Inj-9.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-04-02 15:55:31 : Group 10 of 20, Calling Peaks with MACS2!, 10.93 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C2_Hybrid._.s1UI-10 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C2_Hybrid._.s1UI-10.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-04-02 15:56:27 : Group 11 of 20, Calling Peaks with MACS2!, 11.851 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name hybrid1._.s3UI-11 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/hybrid1._.s3UI-11.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-04-02 15:57:46 : Group 12 of 20, Calling Peaks with MACS2!, 13.165 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name hybrid1._.s3Inj-12 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/hybrid1._.s3Inj-12.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-04-02 15:58:57 : Group 13 of 20, Calling Peaks with MACS2!, 14.347 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name hybrid1._.Other-13 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/hybrid1._.Other-13.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-04-02 16:00:03 : Group 14 of 20, Calling Peaks with MACS2!, 15.457 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name hybrid2._.s3UI-14 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/hybrid2._.s3UI-14.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-04-02 16:01:39 : Group 15 of 20, Calling Peaks with MACS2!, 17.063 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name hybrid2._.1-15 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/hybrid2._.1-15.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-04-02 16:02:21 : Group 16 of 20, Calling Peaks with MACS2!, 17.76 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name resting._.s2UI-16 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/resting._.s2UI-16.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-04-02 16:03:54 : Group 17 of 20, Calling Peaks with MACS2!, 19.302 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name resting._.s1UI-17 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/resting._.s1UI-17.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-04-02 16:05:18 : Group 18 of 20, Calling Peaks with MACS2!, 20.698 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name resting._.s3UI-18 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/resting._.s3UI-18.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-04-02 16:06:21 : Group 19 of 20, Calling Peaks with MACS2!, 21.762 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name resting._.s3Inj-19 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/resting._.s3Inj-19.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-04-02 16:07:06 : Group 20 of 20, Calling Peaks with MACS2!, 22.507 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name resting._.Other-20 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/resting._.Other-20.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## Annotating Peaks : Nearest Gene
## Annotating Peaks : Gene
## Annotating Peaks : TSS
## Annotating Peaks : GC## [1] "/Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/C1_Inj-reproduciblePeaks.gr.rds"## Annotating Peaks : Nearest Gene
## Annotating Peaks : Gene
## Annotating Peaks : TSS
## Annotating Peaks : GC## [1] "/Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/C2_Hybrid-reproduciblePeaks.gr.rds"## Annotating Peaks : Nearest Gene
## Annotating Peaks : Gene
## Annotating Peaks : TSS
## Annotating Peaks : GC## [1] "/Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/hybrid1-reproduciblePeaks.gr.rds"## Annotating Peaks : Nearest Gene
## Annotating Peaks : Gene
## Annotating Peaks : TSS
## Annotating Peaks : GC## [1] "/Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/hybrid2-reproduciblePeaks.gr.rds"## Annotating Peaks : Nearest Gene
## Annotating Peaks : Gene
## Annotating Peaks : TSS
## Annotating Peaks : GC## [1] "/Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/resting-reproduciblePeaks.gr.rds"## Converged after 7 iterations!
## Plotting Ggplot!archProj <- addPeakMatrix(archProj)## ArchR logging to : ArchRLogs/ArchR-addPeakMatrix-13a095ed20471-Date-2021-04-02_Time-16-08-39.log
## If there is an issue, please report to github with logFile!
## 2021-04-02 16:08:40 : Batch Execution w/ safelapply!, 0 mins elapsed.
## .createArrowGroup : Arrow Group already exists! Dropping Group from ArrowFile! This will take ~10-30 seconds!
## .dropGroupsFromArrow : Initializing Temp ArrowFile
## .dropGroupsFromArrow : Adding Metadata to Temp ArrowFile
## .dropGroupsFromArrow : Adding SubGroups to Temp ArrowFile
## .dropGroupsFromArrow : Move Temp ArrowFile to ArrowFile
## 2021-04-02 16:08:47 : Adding s1Inj to PeakMatrix for Chr (1 of 20)!, 0.01 mins elapsed.
## 2021-04-02 16:08:50 : Adding s1Inj to PeakMatrix for Chr (2 of 20)!, 0.053 mins elapsed.
## 2021-04-02 16:08:52 : Adding s1Inj to PeakMatrix for Chr (3 of 20)!, 0.095 mins elapsed.
## 2021-04-02 16:08:55 : Adding s1Inj to PeakMatrix for Chr (4 of 20)!, 0.133 mins elapsed.
## 2021-04-02 16:08:57 : Adding s1Inj to PeakMatrix for Chr (5 of 20)!, 0.176 mins elapsed.
## 2021-04-02 16:09:00 : Adding s1Inj to PeakMatrix for Chr (6 of 20)!, 0.215 mins elapsed.
## 2021-04-02 16:09:02 : Adding s1Inj to PeakMatrix for Chr (7 of 20)!, 0.256 mins elapsed.
## 2021-04-02 16:09:05 : Adding s1Inj to PeakMatrix for Chr (8 of 20)!, 0.297 mins elapsed.
## 2021-04-02 16:09:07 : Adding s1Inj to PeakMatrix for Chr (9 of 20)!, 0.336 mins elapsed.
## 2021-04-02 16:09:09 : Adding s1Inj to PeakMatrix for Chr (10 of 20)!, 0.373 mins elapsed.
## 2021-04-02 16:09:11 : Adding s1Inj to PeakMatrix for Chr (11 of 20)!, 0.411 mins elapsed.
## 2021-04-02 16:09:14 : Adding s1Inj to PeakMatrix for Chr (12 of 20)!, 0.451 mins elapsed.
## 2021-04-02 16:09:16 : Adding s1Inj to PeakMatrix for Chr (13 of 20)!, 0.487 mins elapsed.
## 2021-04-02 16:09:18 : Adding s1Inj to PeakMatrix for Chr (14 of 20)!, 0.523 mins elapsed.
## 2021-04-02 16:09:20 : Adding s1Inj to PeakMatrix for Chr (15 of 20)!, 0.558 mins elapsed.
## 2021-04-02 16:09:22 : Adding s1Inj to PeakMatrix for Chr (16 of 20)!, 0.594 mins elapsed.
## 2021-04-02 16:09:25 : Adding s1Inj to PeakMatrix for Chr (17 of 20)!, 0.629 mins elapsed.
## 2021-04-02 16:09:27 : Adding s1Inj to PeakMatrix for Chr (18 of 20)!, 0.665 mins elapsed.
## 2021-04-02 16:09:29 : Adding s1Inj to PeakMatrix for Chr (19 of 20)!, 0.7 mins elapsed.
## 2021-04-02 16:09:31 : Adding s1Inj to PeakMatrix for Chr (20 of 20)!, 0.734 mins elapsed.
## .createArrowGroup : Arrow Group already exists! Dropping Group from ArrowFile! This will take ~10-30 seconds!
## .dropGroupsFromArrow : Initializing Temp ArrowFile
## .dropGroupsFromArrow : Adding Metadata to Temp ArrowFile
## .dropGroupsFromArrow : Adding SubGroups to Temp ArrowFile
## .dropGroupsFromArrow : Move Temp ArrowFile to ArrowFile
## 2021-04-02 16:09:40 : Adding s3Inj to PeakMatrix for Chr (1 of 20)!, 0.01 mins elapsed.
## 2021-04-02 16:09:43 : Adding s3Inj to PeakMatrix for Chr (2 of 20)!, 0.05 mins elapsed.
## 2021-04-02 16:09:45 : Adding s3Inj to PeakMatrix for Chr (3 of 20)!, 0.089 mins elapsed.
## 2021-04-02 16:09:47 : Adding s3Inj to PeakMatrix for Chr (4 of 20)!, 0.124 mins elapsed.
## 2021-04-02 16:09:50 : Adding s3Inj to PeakMatrix for Chr (5 of 20)!, 0.162 mins elapsed.
## 2021-04-02 16:09:52 : Adding s3Inj to PeakMatrix for Chr (6 of 20)!, 0.2 mins elapsed.
## 2021-04-02 16:09:54 : Adding s3Inj to PeakMatrix for Chr (7 of 20)!, 0.236 mins elapsed.
## 2021-04-02 16:09:56 : Adding s3Inj to PeakMatrix for Chr (8 of 20)!, 0.278 mins elapsed.
## 2021-04-02 16:09:59 : Adding s3Inj to PeakMatrix for Chr (9 of 20)!, 0.324 mins elapsed.
## 2021-04-02 16:10:02 : Adding s3Inj to PeakMatrix for Chr (10 of 20)!, 0.363 mins elapsed.
## 2021-04-02 16:10:04 : Adding s3Inj to PeakMatrix for Chr (11 of 20)!, 0.404 mins elapsed.
## 2021-04-02 16:10:07 : Adding s3Inj to PeakMatrix for Chr (12 of 20)!, 0.445 mins elapsed.
## 2021-04-02 16:10:09 : Adding s3Inj to PeakMatrix for Chr (13 of 20)!, 0.479 mins elapsed.
## 2021-04-02 16:10:11 : Adding s3Inj to PeakMatrix for Chr (14 of 20)!, 0.514 mins elapsed.
## 2021-04-02 16:10:13 : Adding s3Inj to PeakMatrix for Chr (15 of 20)!, 0.551 mins elapsed.
## 2021-04-02 16:10:15 : Adding s3Inj to PeakMatrix for Chr (16 of 20)!, 0.591 mins elapsed.
## 2021-04-02 16:10:17 : Adding s3Inj to PeakMatrix for Chr (17 of 20)!, 0.626 mins elapsed.
## 2021-04-02 16:10:20 : Adding s3Inj to PeakMatrix for Chr (18 of 20)!, 0.663 mins elapsed.
## 2021-04-02 16:10:22 : Adding s3Inj to PeakMatrix for Chr (19 of 20)!, 0.696 mins elapsed.
## 2021-04-02 16:10:24 : Adding s3Inj to PeakMatrix for Chr (20 of 20)!, 0.728 mins elapsed.
## .createArrowGroup : Arrow Group already exists! Dropping Group from ArrowFile! This will take ~10-30 seconds!
## .dropGroupsFromArrow : Initializing Temp ArrowFile
## .dropGroupsFromArrow : Adding Metadata to Temp ArrowFile
## .dropGroupsFromArrow : Adding SubGroups to Temp ArrowFile
## .dropGroupsFromArrow : Move Temp ArrowFile to ArrowFile
## 2021-04-02 16:10:34 : Adding s3UI to PeakMatrix for Chr (1 of 20)!, 0.012 mins elapsed.
## 2021-04-02 16:10:36 : Adding s3UI to PeakMatrix for Chr (2 of 20)!, 0.049 mins elapsed.
## 2021-04-02 16:10:38 : Adding s3UI to PeakMatrix for Chr (3 of 20)!, 0.087 mins elapsed.
## 2021-04-02 16:10:41 : Adding s3UI to PeakMatrix for Chr (4 of 20)!, 0.124 mins elapsed.
## 2021-04-02 16:10:43 : Adding s3UI to PeakMatrix for Chr (5 of 20)!, 0.162 mins elapsed.
## 2021-04-02 16:10:45 : Adding s3UI to PeakMatrix for Chr (6 of 20)!, 0.202 mins elapsed.
## 2021-04-02 16:10:48 : Adding s3UI to PeakMatrix for Chr (7 of 20)!, 0.239 mins elapsed.
## 2021-04-02 16:10:50 : Adding s3UI to PeakMatrix for Chr (8 of 20)!, 0.281 mins elapsed.
## 2021-04-02 16:10:53 : Adding s3UI to PeakMatrix for Chr (9 of 20)!, 0.322 mins elapsed.
## 2021-04-02 16:10:55 : Adding s3UI to PeakMatrix for Chr (10 of 20)!, 0.361 mins elapsed.
## 2021-04-02 16:10:57 : Adding s3UI to PeakMatrix for Chr (11 of 20)!, 0.397 mins elapsed.
## 2021-04-02 16:11:00 : Adding s3UI to PeakMatrix for Chr (12 of 20)!, 0.438 mins elapsed.
## 2021-04-02 16:11:02 : Adding s3UI to PeakMatrix for Chr (13 of 20)!, 0.475 mins elapsed.
## 2021-04-02 16:11:04 : Adding s3UI to PeakMatrix for Chr (14 of 20)!, 0.51 mins elapsed.
## 2021-04-02 16:11:06 : Adding s3UI to PeakMatrix for Chr (15 of 20)!, 0.545 mins elapsed.
## 2021-04-02 16:11:08 : Adding s3UI to PeakMatrix for Chr (16 of 20)!, 0.578 mins elapsed.
## 2021-04-02 16:11:10 : Adding s3UI to PeakMatrix for Chr (17 of 20)!, 0.613 mins elapsed.
## 2021-04-02 16:11:12 : Adding s3UI to PeakMatrix for Chr (18 of 20)!, 0.652 mins elapsed.
## 2021-04-02 16:11:14 : Adding s3UI to PeakMatrix for Chr (19 of 20)!, 0.685 mins elapsed.
## 2021-04-02 16:11:16 : Adding s3UI to PeakMatrix for Chr (20 of 20)!, 0.719 mins elapsed.
## .createArrowGroup : Arrow Group already exists! Dropping Group from ArrowFile! This will take ~10-30 seconds!
## .dropGroupsFromArrow : Initializing Temp ArrowFile
## .dropGroupsFromArrow : Adding Metadata to Temp ArrowFile
## .dropGroupsFromArrow : Adding SubGroups to Temp ArrowFile
## .dropGroupsFromArrow : Move Temp ArrowFile to ArrowFile
## 2021-04-02 16:11:26 : Adding s2UI to PeakMatrix for Chr (1 of 20)!, 0.01 mins elapsed.
## 2021-04-02 16:11:28 : Adding s2UI to PeakMatrix for Chr (2 of 20)!, 0.048 mins elapsed.
## 2021-04-02 16:11:30 : Adding s2UI to PeakMatrix for Chr (3 of 20)!, 0.084 mins elapsed.
## 2021-04-02 16:11:32 : Adding s2UI to PeakMatrix for Chr (4 of 20)!, 0.117 mins elapsed.
## 2021-04-02 16:11:34 : Adding s2UI to PeakMatrix for Chr (5 of 20)!, 0.153 mins elapsed.
## 2021-04-02 16:11:36 : Adding s2UI to PeakMatrix for Chr (6 of 20)!, 0.188 mins elapsed.
## 2021-04-02 16:11:39 : Adding s2UI to PeakMatrix for Chr (7 of 20)!, 0.222 mins elapsed.
## 2021-04-02 16:11:41 : Adding s2UI to PeakMatrix for Chr (8 of 20)!, 0.257 mins elapsed.
## 2021-04-02 16:11:43 : Adding s2UI to PeakMatrix for Chr (9 of 20)!, 0.291 mins elapsed.
## 2021-04-02 16:11:45 : Adding s2UI to PeakMatrix for Chr (10 of 20)!, 0.325 mins elapsed.
## 2021-04-02 16:11:47 : Adding s2UI to PeakMatrix for Chr (11 of 20)!, 0.358 mins elapsed.
## 2021-04-02 16:11:49 : Adding s2UI to PeakMatrix for Chr (12 of 20)!, 0.393 mins elapsed.
## 2021-04-02 16:11:51 : Adding s2UI to PeakMatrix for Chr (13 of 20)!, 0.426 mins elapsed.
## 2021-04-02 16:11:53 : Adding s2UI to PeakMatrix for Chr (14 of 20)!, 0.458 mins elapsed.
## 2021-04-02 16:11:55 : Adding s2UI to PeakMatrix for Chr (15 of 20)!, 0.49 mins elapsed.
## 2021-04-02 16:11:57 : Adding s2UI to PeakMatrix for Chr (16 of 20)!, 0.523 mins elapsed.
## 2021-04-02 16:11:59 : Adding s2UI to PeakMatrix for Chr (17 of 20)!, 0.555 mins elapsed.
## 2021-04-02 16:12:00 : Adding s2UI to PeakMatrix for Chr (18 of 20)!, 0.588 mins elapsed.
## 2021-04-02 16:12:02 : Adding s2UI to PeakMatrix for Chr (19 of 20)!, 0.619 mins elapsed.
## 2021-04-02 16:12:04 : Adding s2UI to PeakMatrix for Chr (20 of 20)!, 0.651 mins elapsed.
## .createArrowGroup : Arrow Group already exists! Dropping Group from ArrowFile! This will take ~10-30 seconds!
## .dropGroupsFromArrow : Initializing Temp ArrowFile
## .dropGroupsFromArrow : Adding Metadata to Temp ArrowFile
## .dropGroupsFromArrow : Adding SubGroups to Temp ArrowFile
## .dropGroupsFromArrow : Move Temp ArrowFile to ArrowFile
## 2021-04-02 16:12:13 : Adding s2Inj to PeakMatrix for Chr (1 of 20)!, 0.01 mins elapsed.
## 2021-04-02 16:12:15 : Adding s2Inj to PeakMatrix for Chr (2 of 20)!, 0.051 mins elapsed.
## 2021-04-02 16:12:18 : Adding s2Inj to PeakMatrix for Chr (3 of 20)!, 0.09 mins elapsed.
## 2021-04-02 16:12:20 : Adding s2Inj to PeakMatrix for Chr (4 of 20)!, 0.126 mins elapsed.
## 2021-04-02 16:12:22 : Adding s2Inj to PeakMatrix for Chr (5 of 20)!, 0.163 mins elapsed.
## 2021-04-02 16:12:24 : Adding s2Inj to PeakMatrix for Chr (6 of 20)!, 0.198 mins elapsed.
## 2021-04-02 16:12:26 : Adding s2Inj to PeakMatrix for Chr (7 of 20)!, 0.232 mins elapsed.
## 2021-04-02 16:12:28 : Adding s2Inj to PeakMatrix for Chr (8 of 20)!, 0.267 mins elapsed.
## 2021-04-02 16:12:30 : Adding s2Inj to PeakMatrix for Chr (9 of 20)!, 0.303 mins elapsed.
## 2021-04-02 16:12:32 : Adding s2Inj to PeakMatrix for Chr (10 of 20)!, 0.337 mins elapsed.
## 2021-04-02 16:12:35 : Adding s2Inj to PeakMatrix for Chr (11 of 20)!, 0.374 mins elapsed.
## 2021-04-02 16:12:37 : Adding s2Inj to PeakMatrix for Chr (12 of 20)!, 0.409 mins elapsed.
## 2021-04-02 16:12:39 : Adding s2Inj to PeakMatrix for Chr (13 of 20)!, 0.442 mins elapsed.
## 2021-04-02 16:12:41 : Adding s2Inj to PeakMatrix for Chr (14 of 20)!, 0.474 mins elapsed.
## 2021-04-02 16:12:43 : Adding s2Inj to PeakMatrix for Chr (15 of 20)!, 0.507 mins elapsed.
## 2021-04-02 16:12:45 : Adding s2Inj to PeakMatrix for Chr (16 of 20)!, 0.539 mins elapsed.
## 2021-04-02 16:12:46 : Adding s2Inj to PeakMatrix for Chr (17 of 20)!, 0.57 mins elapsed.
## 2021-04-02 16:12:48 : Adding s2Inj to PeakMatrix for Chr (18 of 20)!, 0.603 mins elapsed.
## 2021-04-02 16:12:50 : Adding s2Inj to PeakMatrix for Chr (19 of 20)!, 0.634 mins elapsed.
## 2021-04-02 16:12:52 : Adding s2Inj to PeakMatrix for Chr (20 of 20)!, 0.665 mins elapsed.
## .createArrowGroup : Arrow Group already exists! Dropping Group from ArrowFile! This will take ~10-30 seconds!
## .dropGroupsFromArrow : Initializing Temp ArrowFile
## .dropGroupsFromArrow : Adding Metadata to Temp ArrowFile
## .dropGroupsFromArrow : Adding SubGroups to Temp ArrowFile
## .dropGroupsFromArrow : Move Temp ArrowFile to ArrowFile
## 2021-04-02 16:13:01 : Adding s1UI to PeakMatrix for Chr (1 of 20)!, 0.01 mins elapsed.
## 2021-04-02 16:13:03 : Adding s1UI to PeakMatrix for Chr (2 of 20)!, 0.043 mins elapsed.
## 2021-04-02 16:13:05 : Adding s1UI to PeakMatrix for Chr (3 of 20)!, 0.086 mins elapsed.
## 2021-04-02 16:13:08 : Adding s1UI to PeakMatrix for Chr (4 of 20)!, 0.126 mins elapsed.
## 2021-04-02 16:13:10 : Adding s1UI to PeakMatrix for Chr (5 of 20)!, 0.159 mins elapsed.
## 2021-04-02 16:13:12 : Adding s1UI to PeakMatrix for Chr (6 of 20)!, 0.192 mins elapsed.
## 2021-04-02 16:13:14 : Adding s1UI to PeakMatrix for Chr (7 of 20)!, 0.228 mins elapsed.
## 2021-04-02 16:13:16 : Adding s1UI to PeakMatrix for Chr (8 of 20)!, 0.267 mins elapsed.
## 2021-04-02 16:13:18 : Adding s1UI to PeakMatrix for Chr (9 of 20)!, 0.304 mins elapsed.
## 2021-04-02 16:13:20 : Adding s1UI to PeakMatrix for Chr (10 of 20)!, 0.339 mins elapsed.
## 2021-04-02 16:13:22 : Adding s1UI to PeakMatrix for Chr (11 of 20)!, 0.374 mins elapsed.
## 2021-04-02 16:13:25 : Adding s1UI to PeakMatrix for Chr (12 of 20)!, 0.41 mins elapsed.
## 2021-04-02 16:13:27 : Adding s1UI to PeakMatrix for Chr (13 of 20)!, 0.443 mins elapsed.
## 2021-04-02 16:13:29 : Adding s1UI to PeakMatrix for Chr (14 of 20)!, 0.479 mins elapsed.
## 2021-04-02 16:13:31 : Adding s1UI to PeakMatrix for Chr (15 of 20)!, 0.524 mins elapsed.
## 2021-04-02 16:13:34 : Adding s1UI to PeakMatrix for Chr (16 of 20)!, 0.567 mins elapsed.
## 2021-04-02 16:13:36 : Adding s1UI to PeakMatrix for Chr (17 of 20)!, 0.601 mins elapsed.
## 2021-04-02 16:13:38 : Adding s1UI to PeakMatrix for Chr (18 of 20)!, 0.639 mins elapsed.
## 2021-04-02 16:13:41 : Adding s1UI to PeakMatrix for Chr (19 of 20)!, 0.677 mins elapsed.
## 2021-04-02 16:13:43 : Adding s1UI to PeakMatrix for Chr (20 of 20)!, 0.721 mins elapsed.
## Overriding previous entry for ReadsInPeaks
## Overriding previous entry for FRIP
## ArchR logging successful to : ArchRLogs/ArchR-addPeakMatrix-13a095ed20471-Date-2021-04-02_Time-16-08-39.log## get ArchR markers for motif enrichment
markerPeaks_subHybrid <- getMarkerFeatures(
ArchRProj = archProj,
useMatrix = "PeakMatrix",
groupBy = "predAnn",
bias = c("TSSEnrichment", "log10(nFrags)"),
testMethod = "wilcoxon")## ArchR logging to : ArchRLogs/ArchR-getMarkerFeatures-13a0940556392-Date-2021-04-02_Time-16-13-47.log
## If there is an issue, please report to github with logFile!
## MatrixClass = Sparse.Integer.Matrix
## 2021-04-02 16:13:47 : Matching Known Biases, 0.005 mins elapsed.
## 2021-04-02 16:13:51 : Computing Pairwise Tests (1 of 5), 0.069 mins elapsed.
## 2021-04-02 16:14:30 : Computing Pairwise Tests (2 of 5), 0.724 mins elapsed.
## 2021-04-02 16:15:02 : Computing Pairwise Tests (3 of 5), 1.252 mins elapsed.
## 2021-04-02 16:15:30 : Computing Pairwise Tests (4 of 5), 1.721 mins elapsed.
## 2021-04-02 16:15:59 : Computing Pairwise Tests (5 of 5), 2.196 mins elapsed.
## ###########
## 2021-04-02 16:16:27 : Completed Pairwise Tests, 2.676 mins elapsed.
## ###########
## ArchR logging successful to : ArchRLogs/ArchR-getMarkerFeatures-13a0940556392-Date-2021-04-02_Time-16-13-47.logpeakMarkers_subHybrid <- getMarkers(markerPeaks_subHybrid,
cutOff = "FDR <= 0.01 & Log2FC >= 1.25",
returnGR = TRUE)
heatmapMarkerPeaks <- plotMarkerHeatmap(
seMarker = markerPeaks_subHybrid,
cutOff = "FDR <= 0.01 & Log2FC >= 1.25",
transpose = TRUE
)## ArchR logging to : ArchRLogs/ArchR-plotMarkerHeatmap-13a0925fc273d-Date-2021-04-02_Time-16-16-30.log
## If there is an issue, please report to github with logFile!
## Identified 20690 markers!
## Preparing Main Heatmap..
## ArchR logging successful to : ArchRLogs/ArchR-plotMarkerHeatmap-13a0925fc273d-Date-2021-04-02_Time-16-16-30.logComplexHeatmap::draw(heatmapMarkerPeaks, heatmap_legend_side = "bot", annotation_legend_side = "bot")
saveArchRProject(archProj,
outputDirectory = "../data/scATAC/archR_subHybridClusters/")## Copying ArchRProject to new outputDirectory : /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_subHybridClusters
## Copying Arrow Files...
## Copying Arrow Files (1 of 6)
## Copying Arrow Files (2 of 6)
## Copying Arrow Files (3 of 6)
## Copying Arrow Files (4 of 6)
## Copying Arrow Files (5 of 6)
## Copying Arrow Files (6 of 6)
## Getting ImputeWeights
## No imputeWeights found, returning NULL
## Copying Other Files...
## Copying Other Files (1 of 12): Annotations
## Copying Other Files (2 of 12): Embeddings
## Copying Other Files (3 of 12): GroupCoverages
## Copying Other Files (4 of 12): IterativeLSI
## Copying Other Files (5 of 12): PeakCalls
## Copying Other Files (6 of 12): Plots
## Copying Other Files (7 of 12): s1Inj
## Copying Other Files (8 of 12): s1UI
## Copying Other Files (9 of 12): s2Inj
## Copying Other Files (10 of 12): s2UI
## Copying Other Files (11 of 12): s3Inj
## Copying Other Files (12 of 12): s3UI
## Saving ArchRProject...
## Loading ArchRProject...
## Successfully loaded ArchRProject!
##
## / |
## / \
## . / |.
## \\\ / |.
## \\\ / `|.
## \\\ / |.
## \ / |\
## \\#####\ / ||
## ==###########> / ||
## \\##==......\ / ||
## ______ = =|__ /__ || \\\
## ,--' ,----`-,__ ___/' --,-`-===================##========>
## \ ' ##_______ _____ ,--,__,=##,__ ///
## , __== ___,-,__,--'#' ===' `-' | ##,-/
## -,____,---' \\####\\________________,--\\_##,/
## ___ .______ ______ __ __ .______
## / \ | _ \ / || | | | | _ \
## / ^ \ | |_) | | ,----'| |__| | | |_) |
## / /_\ \ | / | | | __ | | /
## / _____ \ | |\ \\___ | `----.| | | | | |\ \\___.
## /__/ \__\ | _| `._____| \______||__| |__| | _| `._____|
##
##
## ___ .______ ______ __ __ .______
## / \ | _ \ / || | | | | _ \
## / ^ \ | |_) | | ,----'| |__| | | |_) |
## / /_\ \ | / | | | __ | | /
## / _____ \ | |\ \\___ | `----.| | | | | |\ \\___.
## /__/ \__\ | _| `._____| \______||__| |__| | _| `._____|
## ## class: ArchRProject
## outputDirectory: /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_subHybridClusters
## samples(6): s1Inj s3Inj ... s2Inj s1UI
## sampleColData names(1): ArrowFiles
## cellColData names(22): Sample TSSEnrichment ... FRIP predAnn
## numberOfCells(1): 4076
## medianTSS(1): 16.2915
## medianFrags(1): 11492.5Color RNA-seq data based on new markers
plotByExpr2 <- function(rna, genes, counts, name="feature", scale=TRUE){
genes <- genes[genes %in% rownames(countHBC)]
if(scale){
scaledCounts <- t(scale(t(counts[genes,]+1e-5)))
scaledCounts <- scaledCounts[!is.na(rowVars(scaledCounts)),]
yhat <- colSums(scaledCounts)
} else {
yhat <- colSums(counts[genes,]+1e-5)
}
rna <- AddMetaData(rna, yhat, col.name = name)
# featureplot
p1 <- FeaturePlot(rna, features = name) + theme(legend.position = "none")
# boxplot
df <- data.frame(y=yhat, cluster=rna$ctHybridSub)
p2 <- ggplot(df, aes(x=cluster, y=yhat)) +
geom_boxplot() +
theme_classic()
cowplot::plot_grid(p1, p2, nrow=1)
}
peakSet <- archProj@peakSet
markersHybrid <- peakMarkers_subHybrid$C2_Hybrid
markersHybrid1 <- peakMarkers_subHybrid$hybrid1
markersHybrid2 <- peakMarkers_subHybrid$hybrid2
peaksHybrid <- peakSet[queryHits(findOverlaps(peakSet, markersHybrid, type = "equal"))]
stopifnot(length(peaksHybrid) == length(markersHybrid))
peaksHybridFunnel <- peakSet[queryHits(findOverlaps(peakSet, markersHybrid1, type = "equal"))]
stopifnot(length(peaksHybridFunnel) == length(markersHybrid1))
peaksHybridHorn <- peakSet[queryHits(findOverlaps(peakSet, markersHybrid2, type = "equal"))]
stopifnot(length(peaksHybridHorn) == length(markersHybrid2))
plotByExpr2(rna,
genes=intersect(unique(mcols(peaksHybrid)$nearestGene), rownames(countHBC)),
counts=countHBC,
name="HybridMarkers",
scale=TRUE)
plotByExpr2(rna,
genes=intersect(unique(mcols(peaksHybridFunnel)$nearestGene), rownames(countHBC)),
counts=countHBC,
name="Hybrid1Markers",
scale=TRUE)
plotByExpr2(rna,
genes=intersect(unique(mcols(peaksHybridHorn)$nearestGene), rownames(countHBC)),
counts=countHBC,
name="Hybrid2Markers",
scale=TRUE)
Motif enrichment for sub-hyrbid clusters
archProj <- addMotifAnnotations(ArchRProj = archProj,
motifSet = "cisbp",
name = "Motif",
force=TRUE)## No methods found in package 'IRanges' for request: 'score' when loading 'TFBSTools'## ArchR logging to : ArchRLogs/ArchR-addMotifAnnotations-13a092ec247f5-Date-2021-04-02_Time-16-17-32.log
## If there is an issue, please report to github with logFile!## peakAnnotation name already exists! Overriding.## 2021-04-02 16:17:35 : Gettting Motif Set, Species : Mus musculus, 0.006 mins elapsed.## Using version 2 motifs!## 2021-04-02 16:17:41 : Finding Motif Positions with motifmatchr!, 0.098 mins elapsed.## 2021-04-02 16:20:59 : Creating Motif Overlap Matrix, 3.392 mins elapsed.## 2021-04-02 16:21:02 : Finished Getting Motif Info!, 3.458 mins elapsed.## ArchR logging successful to : ArchRLogs/ArchR-addMotifAnnotations-13a092ec247f5-Date-2021-04-02_Time-16-17-32.log## upregulated motifs
motifsUp <- peakAnnoEnrichment(
seMarker = markerPeaks_subHybrid,
ArchRProj = archProj,
peakAnnotation = "Motif",
cutOff = "FDR <= 0.05 & Log2FC >= 0.5"
)## ArchR logging to : ArchRLogs/ArchR-peakAnnoEnrichment-13a0955023dce-Date-2021-04-02_Time-16-21-07.log
## If there is an issue, please report to github with logFile!## 2021-04-02 16:21:15 : Computing Enrichments 1 of 5, 0.135 mins elapsed.## 2021-04-02 16:21:16 : Computing Enrichments 2 of 5, 0.151 mins elapsed.## 2021-04-02 16:21:17 : Computing Enrichments 3 of 5, 0.168 mins elapsed.## 2021-04-02 16:21:18 : Computing Enrichments 4 of 5, 0.187 mins elapsed.## 2021-04-02 16:21:20 : Computing Enrichments 5 of 5, 0.209 mins elapsed.## ArchR logging successful to : ArchRLogs/ArchR-peakAnnoEnrichment-13a0955023dce-Date-2021-04-02_Time-16-21-07.log# the top 10 TFs are same for all 3 hybrid clusters which is why there's only 10 cols for all 3 clusters.
heatmapMotifsUp <- plotEnrichHeatmap(motifsUp, transpose = TRUE)## ArchR logging to : ArchRLogs/ArchR-plotEnrichHeatmap-13a09627cafc-Date-2021-04-02_Time-16-21-21.log
## If there is an issue, please report to github with logFile!## Adding Annotations..## Preparing Main Heatmap..ComplexHeatmap::draw(heatmapMotifsUp, heatmap_legend_side = "bot", annotation_legend_side = "bot")
# note that Gm5294 is an ortholog of the FOXL3-OT1 lncRNA and paralog of FOXL1
getTopDf <- function(motifsRes, cl){
df <- data.frame(TF = rownames(motifsRes), mlog10Padj = assay(motifsRes)[,cl])
df <- df[order(df$mlog10Padj, decreasing = TRUE),]
df$rank <- seq_len(nrow(df))
return(df)
}
plotMotifs <- function(motifsRes, cl){
df <- getTopDf(motifsRes, cl)
ggUp <- ggplot(df, aes(rank, mlog10Padj, color = mlog10Padj)) +
geom_point(size = 1) +
ggrepel::geom_label_repel(
data = df[rev(seq_len(30)), ], aes(x = rank, y = mlog10Padj, label = TF),
size = 1.5,
nudge_x = 2,
color = "black"
) + theme_ArchR() +
ylab("-log10(P-adj) Motif Enrichment") +
xlab("Rank Sorted TFs Enriched") +
scale_color_gradientn(colors = paletteContinuous(set = "comet"))
ggUp
}
# the upregulated markers for all hybrid clusters are enriched in Fox motifs.
plotMotifs(motifsUp, "C2_Hybrid") + ggtitle("Hybrid")
plotMotifs(motifsUp, "hybrid1") + ggtitle("hybrid1")
plotMotifs(motifsUp, "hybrid1") + ggtitle("hybrid1")
# Fox motifs are not enriched in either resting or activated cell types.
plotMotifs(motifsUp, "resting") + ggtitle("Resting")
plotMotifs(motifsUp, "C1_Inj") + ggtitle("Activated")
Smarcc1 motifs
Smarcc1 was recently found to regulate Keratin genes and thereby cellular fate in mammalian embryogenesis: https://www.nature.com/articles/s41586-020-2647-4. Here, we find that the markers for the activated cells are most enriched in Smarcc1 motifs. To dig a little deeper, we therefore check which of these marker peaks that may be close to Krt genes have Smarcc1 motifs here.
#detach("package:matrixStats", unload=TRUE)
# get the marker peaks for activated cells
actMarkers <- peakMarkers_subHybrid$C1_Inj
# subset to those that have Smarcc1 motif
mot <- archProj@peakAnnotation$Motif
matches <- readRDS(mot$Matches)
markerHits <- queryHits(findOverlaps(query = SummarizedExperiment::rowRanges(matches),
subject = actMarkers,
type = "equal"))
smarcHits <- which(assays(matches, "matches")[[1]][,grep(x=colnames(assays(matches, "matches")[[1]]), pattern="Smarcc1")])
markerSmarcRanges <- SummarizedExperiment::rowRanges(matches)[intersect(markerHits, smarcHits)]
# check which Krt genes are in there
krtWithMotif <- sort(unname(grep(x=mcols(markerSmarcRanges)$nearestGene, pattern="Krt", value=TRUE)))
krtWithMotif## [1] "Krt13" "Krt14" "Krt16" "Krt18" "Krt19" "Krt20"
## [7] "Krt20" "Krt26" "Krt5" "Krt5" "Krt6a" "Krt84"
## [13] "Krtap11-1" "Krtap2-4" "Krtap4-16"# plot their accessibility and expression
for(kk in 1:length(krtWithMotif)){
gene <- krtWithMotif[kk]
if(gene %in% rownames(rna)){
prna <- FeaturePlot(rna, features=gene)
patac <- FeaturePlot(atac, features=gene) + xlim(c(-6,6))
print(cowplot::plot_grid(prna, patac, nrow=1))
} else {
print(FeaturePlot(atac, features=gene) + xlim(c(-6,6)))
}
}## Warning: Could not find Krt13 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 55 rows containing missing values (geom_point).## Warning: Could not find Krt14 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 55 rows containing missing values (geom_point).
## Warning: Could not find Krt16 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 55 rows containing missing values (geom_point).
## Warning: Could not find Krt18 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 55 rows containing missing values (geom_point).
## Warning: Could not find Krt19 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 55 rows containing missing values (geom_point).
## Warning: Could not find Krt20 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 55 rows containing missing values (geom_point).
## Warning: Could not find Krt20 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 55 rows containing missing values (geom_point).
## Warning: Could not find Krt26 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Warning: Removed 55 rows containing missing values (geom_point).## Warning: Could not find Krt5 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 55 rows containing missing values (geom_point).
## Warning: Could not find Krt5 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 55 rows containing missing values (geom_point).
## Warning: Could not find Krt6a in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 55 rows containing missing values (geom_point).
## Warning: Could not find Krt84 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Warning: Removed 55 rows containing missing values (geom_point).## Warning: Could not find Krtap11-1 in the default search locations, found in
## ACTIVITY assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Warning: Removed 55 rows containing missing values (geom_point).## Warning: Could not find Krtap2-4 in the default search locations, found in
## ACTIVITY assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Warning: Removed 55 rows containing missing values (geom_point).## Warning: Could not find Krtap4-16 in the default search locations, found in
## ACTIVITY assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Warning: Removed 55 rows containing missing values (geom_point).
Check expression of Fox TFs and target genes in scRNA-seq data
Peaks that are more accessible in hybrid cells are upregulated in motifs bound by Fox transcription factors. We expect Fox transcription factors to be higher expressed in hybrid clusters. Since Fox transcription factors generally are repressors, we could expect the target genes of Fox transcription factors to be lower expressed in the hybrid clusters.
Evidence from the literature: - Foxg1 expression is in horizontal basal cells of OE, Foxg1 -/- mice are born with defects in the OE: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2878634/ - Fox family TF motifs are most enriched in olfactory receptors: https://arxiv.org/pdf/1201.2933.pdf - Foxl1 are markers of microvillous cells: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2878634/
# selected Fox TFs as identified with motif enrichment
# top 10 motifs TFs from hybrid cells:
selectedFoxTfs <- c("Floxl1", "Foxa1", "Foxb2", "Foxb1", "Foxc1", "Foxc2", "Foxs1", "Foxd1")
selectedFoxTfs <- intersect(selectedFoxTfs, rownames(countHBC))
countsFoxSelected <- countHBC[selectedFoxTfs,]
foxSumSelected <- colSums(countsFoxSelected)
rna <- AddMetaData(rna, foxSumSelected, col.name = "foxSumSelected")
pFoxSum <- FeaturePlot(rna, features = "foxSumSelected")
df <- data.frame(cl=rna$ctHybridSub,
foxSum=foxSumSelected)
pBoxFoxSum <- ggplot(df, aes(x=cl, y=foxSum)) +
geom_boxplot() +
theme_classic() +
ylab("Sum of selected Fox TFs expression")
cowplot::plot_grid(pFoxSum, pBoxFoxSum)
#ggsave("../plots/foxSumRNA.png", width=8)
## scaled
scaledCountsFoxSelected <- t(scale(t(countHBC[selectedFoxTfs,]+1e-5)))
scaledCountsFoxSelected <- scaledCountsFoxSelected[!is.na(rowVars(scaledCountsFoxSelected)),]
foxSumScaledSelected <- colSums(scaledCountsFoxSelected)
rna <- AddMetaData(rna, foxSumScaledSelected, col.name = "foxSumSelected_scaled")
pFoxSumScaled <- FeaturePlot(rna, features = "foxSumSelected_scaled")
df <- data.frame(cl=rna$ctHybridSub,
foxScaled=foxSumScaledSelected)
pBoxFoxSumScaled <- ggplot(df, aes(x=cl, y=foxScaled)) +
geom_boxplot() +
theme_classic() +
ylab("Sum of selected Fox TFs scaled expression")
cowplot::plot_grid(pFoxSumScaled, pBoxFoxSumScaled)
# ggsave("../plots/foxSumRNAScaled.png", width=8)Composite plot
## DR of RNA
pRNA1_fin <- pRNA1 + ggtitle("scRNA-seq data") +
theme(legend.title = element_text(size = 12),
legend.text = element_text(size = 10),
legend.key.size = unit(1/2, "cm"))
## DR of RNA with ATAC markers expression highlighted
pMarkersATAC <- cowplot::plot_grid(plotlist=list(pc1,pc2,pc3), nrow=1)
pMarkersATAC
## RNA and ATAC cell type transfer
pCTrans <- cowplot::plot_grid(ppred4_2 + NoLegend(), ppred2_2, nrow=1)
## markers and motifs for subhybrid clusters
## marker
heatmapMarkerPeaksMatrix <- plotMarkerHeatmap(
seMarker = markerPeaks_subHybrid,
cutOff = "FDR <= 0.01 & Log2FC >= 1.25",
transpose = TRUE,
returnMat = TRUE
)## ArchR logging to : ArchRLogs/ArchR-plotMarkerHeatmap-13a095fd2bfdc-Date-2021-04-02_Time-16-22-06.log
## If there is an issue, please report to github with logFile!## Identified 20690 markers!## ArchR logging successful to : ArchRLogs/ArchR-plotMarkerHeatmap-13a095fd2bfdc-Date-2021-04-02_Time-16-22-06.logpMarker <- Heatmap(heatmapMarkerPeaksMatrix,
cluster_rows = FALSE,
cluster_columns = FALSE,
col = paletteContinuous(set = "solarExtra", n = 100),
show_column_names = FALSE,
row_names_gp = gpar(fontsize = 13),
heatmap_legend_param = list(title = "Z-score"))
## motif
heatmapMotifsMatrix <- plotEnrichHeatmap(motifsUp, transpose = TRUE, returnMatrix = TRUE)## ArchR logging to : ArchRLogs/ArchR-plotEnrichHeatmap-13a0953235aa-Date-2021-04-02_Time-16-22-08.log
## If there is an issue, please report to github with logFile!colnames(heatmapMotifsMatrix) <- unlist(lapply(strsplit(colnames(heatmapMotifsMatrix), split="_"), "[[", 1))
pMotif <- Heatmap(heatmapMotifsMatrix,
cluster_rows = FALSE,
cluster_columns = FALSE,
col = paletteContinuous(set = "comet", n = 100),
column_names_gp = gpar(fontsize = 11),
row_names_gp = gpar(fontsize = 13),
heatmap_legend_param = list(title = "-log10\n p-adj."))
pHeat <- cowplot::plot_grid(grid::grid.grabExpr(ComplexHeatmap::draw(pMarker)),
grid::grid.grabExpr(ComplexHeatmap::draw(pMotif)),
nrow = 2)
pComp1 <- cowplot::plot_grid(pRNA1_fin, pMarkersATAC, nrow=1,
rel_widths = c(0.4, 1), labels = c("a", "b"))
pComp2 <- cowplot::plot_grid(pCTrans, pHeat, ncol = 1,
rel_heights = c(0.8, 1), labels = c("c", "d"))
cowplot::plot_grid(pComp1, pComp2, nrow = 2,
rel_widths= c(1/3, 1))
# ggsave("../plots/integrationComposite.pdf", height=13, width=12)
# ggsave("../plots/integrationComposite.png", height=13, width=12)old plots
Explore the activated HBC RNA-seq clusters
Macrophage markers
FeaturePlot(atac, features = "Msr1")## Warning: Could not find Msr1 in the default search locations, found in ACTIVITY
## assay instead
FeaturePlot(rna, features = "Msr1")
Cell cycle genes
download.file("https://raw.githubusercontent.com/rufletch/p63-HBC-diff/master/ref/cellCycle40.txt", destfile = "~/tmp/ccGenes40.txt")
cc40Genes <- read.table("~/tmp/ccGenes40.txt", stringsAsFactors = FALSE)[,1]
countsCC40 <- countHBC[cc40Genes,]+1e-5
cc40Sum <- colSums(countsCC40)
rna <- AddMetaData(rna, cc40Sum, col.name = "cc40Sum")
FeaturePlot(rna, features = "cc40Sum")
## scaled
scaledCountsCC40 <- t(scale(t(countHBC[cc40Genes,]+1e-5)))
scaledCountsCC40 <- scaledCountsCC40[!is.na(rowVars(scaledCountsCC40)),]
cc40Sum <- colSums(scaledCountsCC40)
rna <- AddMetaData(rna, cc40Sum, col.name = "cc40Sum_scaled")
FeaturePlot(rna, features = "cc40Sum_scaled")
download.file("https://raw.githubusercontent.com/rufletch/p63-HBC-diff/master/ref/cell_cycle.txt", destfile = "~/tmp/ccGenes.txt")
ccGenes <- read.table("~/tmp/ccGenes.txt", stringsAsFactors = FALSE)[,1]
ccGenes <- ccGenes[ccGenes %in% rownames(countHBC)]
countsCC <- countHBC[ccGenes,]+1e-5
ccSum <- colSums(countsCC)
rna <- AddMetaData(rna, ccSum, col.name = "ccSum")
FeaturePlot(rna, features = "ccSum")
## scaled
scaledCountsCC <- t(scale(t(countHBC[ccGenes,]+1e-5)))
scaledCountsCC <- scaledCountsCC[!is.na(rowVars(scaledCountsCC)),]
ccSum <- colSums(scaledCountsCC)
rna <- AddMetaData(rna, ccSum, col.name = "ccSum_scaled")
FeaturePlot(rna, features = "ccSum_scaled")
Wound response genes
Genes obtained from paper on HBC injury, Fig 3D caption: `Krt6a, Krt16, Sprr1a, Sprr2a3, Krtdap, Dmkn, Sbsn, and Hbegf’
wrGenes <- c("Krt6a", "Krt16", "Sprr1a", "Sprr2a3", "Krtdap", "Dmkn", "Sbsn", "Hbegf")
countsWR <- countHBC[wrGenes,]+1e-5
wrSum <- colSums(countsWR)
rna <- AddMetaData(rna, wrSum, col.name = "wrSum")
FeaturePlot(rna, features = "wrSum")
## scaled
scaledCountsWR <- t(scale(t(countHBC[wrGenes,]+1e-5)))
scaledCountsWR <- scaledCountsWR[!is.na(rowVars(scaledCountsWR)),]
wrSum <- colSums(scaledCountsWR)
rna <- AddMetaData(rna, wrSum, col.name = "wrSum_scaled")
FeaturePlot(rna, features = "wrSum_scaled")
Respiratory markers
respGenesRebecca <- c("Gp2", "Cyp4a12b", "Kcnj16", "Tff2", "Chad", "Acsm1", "Kcne3", "Vtcn1", "Kcnj15", "Cldn8", "Pigr", "Muc5b", "Calml3", "Kng2", "Cyp2b10")
## note that Cxcl14 and Meg3 are only described for human samples in the Methods.
respGenesDavid <- c("Foxj1", "Muc5ac", "Cxcl14", "Meg3")
## genes used by Boying
respGenesBoying <- c("Reg3g", "Muc5ac", "Muc5b", "Meg3", "Krt5", "Krt14", "Trp63")
FeaturePlot(rna, features = respGenesRebecca[1:9])
FeaturePlot(rna, features = respGenesRebecca[10:15])
FeaturePlot(rna, features = respGenesDavid)
FeaturePlot(rna, features = respGenesBoying) + NoLegend()
Session Info
sessionInfo()## R version 3.6.2 (2019-12-12)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.5
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] grid parallel stats4 stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] chromVARmotifs_0.2.0 motifmatchr_1.8.0
## [3] circlize_0.4.8 ComplexHeatmap_2.3.1
## [5] presto_1.0.0 nabor_0.5.0
## [7] BSgenome.Mmusculus.UCSC.mm10_1.4.0 BSgenome_1.54.0
## [9] rtracklayer_1.46.0 Biostrings_2.54.0
## [11] XVector_0.26.0 liger_0.5.0
## [13] ggalluvial_0.12.0 ArchR_0.9.5
## [15] magrittr_1.5 rhdf5_2.30.1
## [17] data.table_1.12.8 patchwork_1.0.0
## [19] Seurat_3.1.5 Signac_0.2.5
## [21] harmony_1.0 Rcpp_1.0.6
## [23] scater_1.14.6 uwot_0.1.5
## [25] scran_1.14.5 forcats_0.4.0
## [27] stringr_1.4.0 dplyr_1.0.3
## [29] purrr_0.3.3 readr_1.3.1
## [31] tidyr_1.0.2 tibble_2.1.3
## [33] tidyverse_1.3.0 irlba_2.3.3
## [35] Matrix_1.3-2 gridExtra_2.3
## [37] wesanderson_0.3.6 pheatmap_1.0.12
## [39] ggplot2_3.3.2 RColorBrewer_1.1-2
## [41] clusterExperiment_2.6.1 bigmemory_4.5.36
## [43] rgl_0.100.30 cowplot_1.0.0
## [45] SingleCellExperiment_1.8.0 SummarizedExperiment_1.16.1
## [47] DelayedArray_0.12.2 BiocParallel_1.20.1
## [49] matrixStats_0.56.0 Biobase_2.46.0
## [51] GenomicRanges_1.38.0 GenomeInfoDb_1.22.0
## [53] IRanges_2.20.2 S4Vectors_0.24.3
## [55] BiocGenerics_0.32.0 slingshot_1.4.0
## [57] princurve_2.1.4
##
## loaded via a namespace (and not attached):
## [1] rsvd_1.0.2 ica_1.0-2
## [3] zinbwave_1.11.6 Rsamtools_2.2.1
## [5] foreach_1.4.7 lmtest_0.9-37
## [7] rprojroot_1.3-2 crayon_1.3.4
## [9] MASS_7.3-51.4 nlme_3.1-142
## [11] backports_1.1.5 reprex_0.3.0
## [13] rlang_0.4.10 ROCR_1.0-7
## [15] readxl_1.3.1 limma_3.42.1
## [17] phylobase_0.8.6 rjson_0.2.20
## [19] CNEr_1.20.0 manipulateWidget_0.10.0
## [21] bit64_0.9-7 glue_1.4.2
## [23] poweRlaw_0.70.2 rngtools_1.5
## [25] sctransform_0.3.2 vipor_0.4.5
## [27] AnnotationDbi_1.48.0 VGAM_1.1-2
## [29] haven_2.2.0 tidyselect_1.1.0
## [31] fitdistrplus_1.0-14 XML_3.99-0.3
## [33] zoo_1.8-7 GenomicAlignments_1.22.1
## [35] xtable_1.8-4 evaluate_0.14
## [37] bibtex_0.4.2.2 cli_2.0.1
## [39] zlibbioc_1.32.0 rstudioapi_0.11
## [41] miniUI_0.1.1.1 locfdr_1.1-8
## [43] shiny_1.4.0 BiocSingular_1.2.1
## [45] xfun_0.12 clue_0.3-57
## [47] cluster_2.1.0 caTools_1.18.0
## [49] doSNOW_1.0.18 ggfittext_0.9.0
## [51] KEGGREST_1.26.1 ggrepel_0.8.1
## [53] ape_5.3 listenv_0.8.0
## [55] TFMPvalue_0.0.8 png_0.1-7
## [57] future_1.16.0 withr_2.1.2
## [59] bitops_1.0-6 plyr_1.8.5
## [61] cellranger_1.1.0 dqrng_0.2.1
## [63] pillar_1.4.3 RcppParallel_4.4.4
## [65] gplots_3.0.1.2 GlobalOptions_0.1.1
## [67] fs_1.3.1 kernlab_0.9-29
## [69] hdf5r_1.3.1 GetoptLong_0.1.8
## [71] DelayedMatrixStats_1.8.0 vctrs_0.3.6
## [73] ellipsis_0.3.0 generics_0.0.2
## [75] NMF_0.21.0 tools_3.6.2
## [77] rncl_0.8.3 beeswarm_0.2.3
## [79] munsell_0.5.0 fastmap_1.0.1
## [81] compiler_3.6.2 httpuv_1.5.2
## [83] pkgmaker_0.31 plotly_4.9.1
## [85] GenomeInfoDbData_1.2.2 edgeR_3.28.0
## [87] riverplot_0.6 lattice_0.20-38
## [89] snow_0.4-3 later_1.0.0
## [91] jsonlite_1.6.1 scales_1.1.0
## [93] pbapply_1.4-2 genefilter_1.68.0
## [95] lazyeval_0.2.2 promises_1.1.0
## [97] doParallel_1.0.15 R.utils_2.9.2
## [99] reticulate_1.14 rmarkdown_2.1
## [101] statmod_1.4.33 webshot_0.5.2
## [103] Rtsne_0.15 softImpute_1.4
## [105] igraph_1.2.4.2 HDF5Array_1.14.1
## [107] survival_3.1-8 yaml_2.2.1
## [109] htmltools_0.5.1.1 memoise_1.1.0
## [111] locfit_1.5-9.1 here_0.1
## [113] viridisLite_0.3.0 digest_0.6.27
## [115] assertthat_0.2.1 mime_0.9
## [117] registry_0.5-1 npsurv_0.4-0
## [119] bigmemory.sri_0.1.3 RSQLite_2.2.0
## [121] future.apply_1.4.0 lsei_1.2-0
## [123] blob_1.2.1 R.oo_1.23.0
## [125] RNeXML_2.4.2 TFBSTools_1.22.0
## [127] splines_3.6.2 labeling_0.3
## [129] Rhdf5lib_1.8.0 Cairo_1.5-10
## [131] RCurl_1.98-1.1 broom_0.5.4
## [133] hms_0.5.3 modelr_0.1.5
## [135] colorspace_1.4-1 ggbeeswarm_0.6.0
## [137] shape_1.4.4 mclust_5.4.5
## [139] RANN_2.6.1 ggseqlogo_0.1
## [141] fansi_0.4.1 R6_2.4.1
## [143] ggridges_0.5.2 lifecycle_0.2.0
## [145] gdata_2.18.0 leiden_0.3.4
## [147] howmany_0.3-1 RcppAnnoy_0.0.16
## [149] iterators_1.0.12 htmlwidgets_1.5.1
## [151] crosstalk_1.0.0 seqLogo_1.50.0
## [153] rvest_0.3.5 globals_0.12.5
## [155] codetools_0.2-16 lubridate_1.7.4
## [157] GO.db_3.10.0 FNN_1.1.3
## [159] gtools_3.8.1 prettyunits_1.1.1
## [161] dbplyr_1.4.2 gridBase_0.4-7
## [163] RSpectra_0.16-0 R.methodsS3_1.7.1
## [165] gtable_0.3.0 tsne_0.1-3
## [167] DBI_1.1.0 httr_1.4.1
## [169] KernSmooth_2.23-16 stringi_1.4.5
## [171] progress_1.2.2 reshape2_1.4.3
## [173] farver_2.0.3 uuid_0.1-2
## [175] annotate_1.64.0 viridis_0.5.1
## [177] gggenes_0.4.0 xml2_1.3.2
## [179] BiocNeighbors_1.4.1 ade4_1.7-13
## [181] bit_1.1-15.2 pkgconfig_2.0.3
## [183] DirichletMultinomial_1.26.0 knitr_1.28